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ABSTRACT 

With  the  U.S.  Navy’s  eontinued  foeus  on  Integrated  Eight  Thru  Power  (lETP)  there  has  been  an 
ever  increasing  effort  to  ensure  an  electrical  distribution  system  that  maintains  maximum 
capabilities  in  the  event  of  system  faults.  This  is  to  ensure  that  the  crew  has  the  ability  to 
complete  real  time  tactical  missions  in  the  event  of  battle  damage  to  any  localized  portions  of  the 
electrical  distribution  system.  Eault  isolation  is  a  priority  component  of  the  U.S.  Navy’s  Next 
Generation  Integrated  Power  System  (NGIPS)  Roadmap,  which  lays  out  the  framework  as  well 
as  milestone  dates  for  future  development.  Non-Intrusive  Toad  Monitoring  (NIEM),  which  has 
been  used  extensively  for  condition  based  maintenance  applications,  could  simultaneously  be 
used  to  enhance  the  existing  zonal  protection  system  employed  with  Multi-Eunction  Monitors 
(MEM).  NIEM  may  be  able  to,  inexpensively,  use  the  existing  current  and  voltage  sensors 
available  from  the  MEM  hardware  to  determine  electrical  loading  which  could  allow  for  faster 
fault  isolation  capability. 

A  test  platform  with  three  5000  watt  synchronous  generators  is  being  constructed  to  emulate  a 
U.S.  Navy  DDG  51  EET  IIA  class  ship  electric  plant.  This  is  being  accomplished  in  order  to 
evaluate  the  feasibility  of  improving  the  fault  isolation  capabilities  of  the  MEM  with  NIEM 
implementation.  The  first  step  in  this  endeavor  will  be  to  electrically  relate  the  test  platform  to 
the  DDG  electric  plant.  In  order  to  accomplish  this  step,  the  fault  simulation  results  from  the  test 
platform  will  be  compared  to  simulated  faults  using  U.S.  Navy  data  from  DDG  51  electric  plants. 
This  will  allow  for  the  fault  isolation  results  from  the  test  platform  to  be  related  to  the  DDG  5 1 
electric  plant. 


Thesis  Supervisor:  Steven  B.  Eeeb 

Title:  Professor  of  Electrical  Engineering  and  Computer  Science  &  Mechanical  Engineering 
Thesis  Supervisor:  Robert  W.  Cox 

Title:  Assistant  Professor  of  Electrical  and  Computer  Engineering,  UNC  Charlotte 


3 


Acknowledgements 


The  author  would  like  to  aeknowledge  the  following  individuals  for  their  assistanee.  Without 
them  this  thesis  would  not  have  been  possible. 

•  Professor  Steven  Leeb  for  providing  me  with  a  ehallenging  and  rewarding  experienee 

•  Professor  Robert  Cox  for  his  enthusiasm,  assistanee  and  feedbaek 

•  Ward  Wiehakool  for  his  teehnioal  expertise  and  willingness  to  help 

•  Uzoma  Orji  for  his  software  and  teehnieal  expertise  and  willingness  to  help 

•  Captain  Norbert  Doerry  for  providing  all  the  information  and  eontaets  neeessary  to  help  with 
my  thesis 

•  Timothy  Zavertnik  for  his  willingness  to  provide  data  and  answer  questions 

•  Keith  Douglas  and  Greg  Elkins  for  their  teehnieal  aeumen  and  ability  to  keep  things  fun 
during  shipboard  installations. 


Finally,  thanks  and  love  to  my  wife  and  daughter  who  have  been  patient  and  supportive 
throughout  this  proeess. 


4 


Table  of  Contents 


Acknowledgements . 4 

Table  of  Contents . 5 

List  of  Figures . 7 

List  of  Tables . 9 

1.0  Introduction . 10 

1.1  Background . 10 

1.2  Non-Intrusive  Load  Monitor  (NILM) . 12 

1.3  Thesis  Objeetives . 13 

2.0  Zonal  Eleetrical  Distribution  System  with  MFM  III . 15 

2.1  High  Speed  Relay  (HSR)  Algorithm . 17 

3.0  Synchronous  Machines . 20 

3.1  Basie  Maehine  Layout . 20 

3.2  Direct  and  Quadrature  Transformation . 21 

3.3  Equal  Mutuals  Model  Description . 25 

4.0  Preliminary  Data  Colleetion . 28 

4.1  DC  test . 28 

4.2  Open-Circuit  Test:  Magnetization  Characteristic  and  Eield  Resistance . 29 

4.3  Short  Circuit  Characteristics . 31 

4.4  Short  Circuit  Testing . 34 

4.4.1  Test  Setup . 35 

4.4.2  Test  Results . 37 

4.4.3  Analysis  of  Data . 40 

5.0  Generator  Modeling . 43 


5 


5.1  Equal  Mutuals  Modeling . 43 

5.1.1  Bench  Generator  Model . 46 

5.1.2  Sensitivity  Analysis . 52 

5.1.2  DDG  Model . 54 

5.2  DDG  and  Bench  Generator  Comparison . 57 

5.2.1  Armature  Reaction  for  DDG  and  Bench  Generator . 63 

6.0  Future  Studies . 65 

6.1  Additional  Work . 65 

6.2  Conclusion . 66 

BIBLIOGRAPHY . 68 

APPENDIX  1  -  MATLAB  Evaluation  of  Xd . 69 

APPENDIX  2  -  MATLAB  tachometer  to  frequency  code . 73 

APPENDIX  3  -  MATLAB  Bench  Generator  Model  Code . 75 

APPENDIX  4  -  MATLAB  DDG-51  Simulation  Code  and  DDG-II  Results . 81 

APPENDIX  5  -  SIMULATION  USERS  GUIDE . 87 

APPENDIX  6  -  MATLAB  MODEL  VERIFICATION . 89 

APPENDIX  7  -  CALCULATION  OF  NAVSEA  DDG  PARAMETERS . 100 


6 


List  of  Figures 

Figure  1:  NILM  architecture  (3) . 12 

Figure  2:  NILM  installation  (3) . 13 

Figure  3:  Ming  Dong  5000  Watt  Bench  Generator . 14 

Figure  4:  Multi  Function  Module  (MFM  III) . 16 

Figure  5:  Electric  Plant  Layout  (DDG51  FLT  IIA) . 17 

Figure  6:  Four  pole,  Y-Connected  Synchronous  Generator  (5) . 20 

Figure  7:  Direct  Axis  Model . 26 

Figure  8:  Quadrature  Axis  Model . 26 

Figure  9;  Bench  Generator  Y-configuration(7) . 29 

Figure  10;  Bench  Generator  OCC . 30 

Figure  11:  SCC  and  OCC  test  setup  (7) . 31 

Figure  12;  Bench  Generator  Iscc . 32 

Figure  13;  Bench  Generator  OCC  and  SCC . 32 

Figure  14;  Bench  Generator  D-Axis  Reactance . 33 

Figure  15:  Short  circuit  envelope  (classical  equation)  (8) . 35 

Figure  16:  Bench  Generator  Prime  Mover  (1.5  HP) . 36 

Figure  17:  Switching  Apparatus  for  short  circuit  implementation . 36 

Figure  18;  Short  Circuit  Results  from  Oscilloscope . 37 

Figure  19;  Field  Current  during  short  circuit . 38 

Figure  20:  Speed  profile  for  Short  Circuit  (Tachometer  and  Inertia  Constant) . 39 

Figure  21;  Evaluation  of  Xd' . 40 

Eigure  22:  Xd"  Determination . 41 

Eigure  23:  Direct  and  Quadrature  Axes  (6) . 45 

Eigure  24;  Bench  Generator  Simulation . 5 1 

Eigure  25;  Bench  Generator  (Simulation,  Classical,  and  Actual)  Results . 52 

Eigure  26:  DDG-I  Short  Circuit  (Initially  Unloaded) . 56 

Eigure  27;  DDG-I  Short  Circuit  (Initially  Eully  Loaded) . 56 

Eigure  28:  Pole  Zero  Plots  for  Bench  Generator  and  DDG-I  Model . 58 

Eigure  29;  Impedance  Bode  Plots  (Bench  Generator  and  DDG-I) . 59 


7 


Figure  30;  DDG-51  Current  Threshold  after  Three  Phase  Fault  (Initially  Unloaded  -  Left), 

(Initially  Fully  Loaded  -  Right) . 60 

Figure  31:  Beneh  Generator  Current  Threshold  (Initially  Unloaded  -  Left),  (Initially  Fully 

Loaded  -  Right) . 60 

Figure  32:  Bench  Generator  Actual  Phase  Current  Crossing  150%  Threshold . 62 

Figure  33:  Effect  of  Armature  Resistance  on  Short  Circuit  Current  for  DDG . 64 

Figure  34:  DDG-II  Short  Circuit  Simulation  (Initially  Unloaded) . 85 

Figure  35:  DDG-II  Short  Circuit  Simulation  Frequency  Response  (Initially  Unloaded) . 86 

Figure  36;  Krause  Simulation  used  to  verify  Equal  Mutuals  Simulation  Model . 89 

Eigure  37:  Krause  Simulation  Results . 90 

Eigure  38;  Equal  Mutuals  Model  Results  for  Krause  Simulation  (Phase,  Direct,  and  Quadrature 

Currents) . 97 

Eigure  39;  Equal  Mutuals  Model  Results  for  Krause  Simulation  (Erequency  and  Degrees) . 98 

Eigure  40;  Equal  Mutuals  Model  Results  for  Krause  Simulation  (Eield  Current  and  Torque)  ....  99 


8 


List  of  Tables 


Table  1;  Bench  Generator  Specifications . 28 

Table  2:  Bench  Generator  Parameters . 48 

Table  3:  Simulation  and  Classical  results  compared  to  Actual  Data . 52 

Table  4:  Bench  Generator  Model  Sensitivity  Parameters  (First  Peak  Current) . 53 

Table  5:  Bench  Generator  Sensitivity  Parameters  (Decay  Envelope) . 54 

Table  6;  DDG-51  generator  parameters  provided  by  NAVSEA . 55 

Table  7:  Generator  Times  for  Eault  Detection  (150%  Current) . 62 

Table  8;  Armature  time  constants  for  DDG  and  Bench  Generator . 63 


9 


Chapter  1 


1.0  Introduction 

The  DDG-51  FLT  IIA  class  ships  currently  use  Multi-function  Monitors  (MFMs)  to  quickly 
locate  and  isolate  electrical  faults.  These  systems  perform  fault  detection  using  the  high  speed 
relay  (HSR)  algorithm  that  monitors  the  magnitude  and  phase  of  the  various  system  voltages  (1). 
One  difficulty  encountered  when  implementing  the  HSR  algorithm  is  the  need  to  be  able  to 
distinguish  between  true  faults  and  momentary  voltage  “spikes”  that  occur  because  of  switching 
events  and  other  normal  activities.  Because  most  spikes  dampen  out  quickly,  they  can  be  easily 
detected  using  additional  samples  or  removed  using  low-pass  fdters.  Both  of  these  options  add 
delays  to  the  HSR  algorithm,  thus  limiting  its  ability  to  detect  true  faults.  Careful  adjustment  of 
the  algorithm’s  thresholds  is  necessary  both  to  prevent  voltage  degradations  that  can  disrupt 
power  quality  during  faults  and  to  avoid  false  alarms.  As  NAVSEA  has  found,  it  is  difficult  to 
determine  optimal  threshold  settings  using  only  laboratory  test  circuits.  Instead,  appropriate 
thresholds  can  only  be  determined  if  one  knows  the  maximum  allowable  waveform  distortion, 
and  this  can  be  different  depending  on  the  shipboard  environment  (2).  One  way  to  solve  this 
problem  is  to  integrate  the  MFM  with  a  tool  known  as  the  non-intrusive  load  monitor  (NILM). 
The  NILM  is  a  device  that  can  classify  various  transient  events  in  power  systems,  and  it  has  been 
used  repeatedly  for  load-identification  in  shipboard  environments  (3).  The  NILM  may  be  able  to 
improve  the  HSR  performance  by  identifying  switching  events  and  by  maintaining  a  list  of 
currently  operating  loads.  With  this  information,  it  becomes  possible  to  dynamically  select 
settings  based  on  current  conditions. 

1.1  Background 

With  the  US  Navy  warships  evolution  to  the  Integrated  Power  Systems  (IPS)  architecture  an 
emphasis  has  been  placed  on  Zonal  Electrical  Distribution  Systems  (ZEDS).  ZEDS  focus  on 
classifying  electrical  loads  into  three  categories  (uninterruptible,  short  term  interrupt,  and  long 
term  interrupt)  (4).  Uninterruptible  loads  would  include  items  such  as  that  necessary  to  execute 
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the  ships  mission  (including  weaponry,  sensors,  communication,  and  propulsion  systems).  IPS 
will  provide  electric  propulsion  and  ship  service  electrical  loads,  from  a  common  power  source, 
into  one  distribution  system.  Historically  this  has  not  been  the  case.  Propulsion  requirements 
have  been  achieved  mechanically  through  either  a  diesel  or  gas  turbine  coupled  to  a  shaft 
through  an  intermediary  reduction  gear.  This  in  turn  has  buffered  the  propulsion  plant  from 
having  any  direct  degradation  effects  as  a  result  of  most  electric  plant  casualties.  Electric  plants 
have  been  maintained  in  redundant  split  bus  configurations  as  well  as  having  separate  vital  and 
non-vital  busses.  The  vital  busses  would  contain  items  necessary  for  war  fighting  as  well  as 
equipment  necessary  for  maintaining  propulsion  such  as  propulsion  lube  oil.  In  the  event  of 
casualties,  the  plant  has  been  designed  to  maximize  continuity  of  power  to  the  vital  bus  while 
also  ensuring  that  the  non-vital  bus  is  the  first  to  be  disconnected.  This  has  been  accomplished 
by  having  breakers  that  open  on  such  items  as  over-current,  under-frequency,  over-frequency, 
under-voltage,  along  with  other  very  specific  criteria.  Electric  plants  have  not  been  designed  to 
logically  assess  where  a  fault  is  coming  from  and  to  only  isolate  power  locally  to  those  affected 
areas.  However,  with  the  advent  of  electric  propulsion  it  is  now  necessary  to  be  able  to  localize 
and  isolate  faults  to  only  those  zones  affected  in  order  to  ensure  maximized  propulsion 
capabilities  as  well  as  war  fighting  capabilities.  This  is  known  as  Integrated  Eight  Through 
Power  (lETP)  and  is  accomplished  by  having  a  ZEDS  that  uses  Multi-Eunction  Monitors  (MEM) 
and  a  High  Speed  Relay  (HSR)  algorithm.  The  MEM  is  used  to  identify  where  faults  are 
originating  and  to  only  isolate  those  areas.  Currently  MEM  III  operation  requires  a  fault 
threshold  to  be  set  which  will  discriminate  faults  from  typical  plant  transients.  The  plant 
transients  must  be  assessed  for  the  worst  case  loading  scenarios  which  could  delay  true  fault 
responses.  Information  obtained  from  NIEM  could  potentially  be  used  to  optimize  fault 
thresholds  based  on  current  plant  loading. 

The  Next  Generation  Integrated  Power  System  (NGIPS)  has  listed  as  challenges  for  future 
implementation  of  AC  ZEDS  the  implementation  of  mission  priority  load  shedding.  NIEM 
could  aid  this  process  by  having  knowledge  of  current  power  loading  as  well  as  the  capability  of 
knowing  which  loads  are  currently  operating  real  time.  Through  this  knowledge  NIEM  could 
communicate  the  need  to  turn  off  non-vital  loads  when  power  is  mission  limiting. 
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The  NGIPS  description  of  what  is  desired  in  Power  Control  Modules  (PCON)  leads  to  a 
conclusion  that  NILM  could  be  developed  to  meet  those  requirements.  The  primary  functions  of 
PCON  which  could  be  improved  by  NILM  are  the  following  (4): 

1 .  Remote  monitoring  and  control  of  NGIPS  modules. 

2.  Mission  priority  load  shedding. 

3.  Fault  detection  and  isolation. 

4.  Maintenance  support. 

1.2  Non-Intrusive  Load  Monitor  (NILM) 

The  NILM  is  a  device  that  measures  at  least  one  voltage,  and  one  line  current  with  a  transducer 
at  a  particular  point  in  an  electrical  distribution  system.  These  measurements  are  used  to 
calculate  both  real  and  reactive  power  envelopes  as  well  as  harmonic  frequency  content.  These 
measurements  and  calculated  values  can  be  used  for  monitoring  which  loads  are  currently 
operating  and  in  what  mode  (e.g.  slow  speed  or  fast  speed  pumps),  as  well  as  for  determining  if 
equipment  is  operating  correctly  or  is  need  of  maintenance.  NILM  can  be  used  for  power 
monitoring  and  evaluation  of  fault  conditions  in  conjunction  with  MFM  and  HSR  algorithms. 


Command  Outputs  Status  Reports 


Operator 

Interface 


Figure  1:  NILM  architecture  (3) 
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The  NILM  architecture  is  shown  in  Figure  1 .  It  allows  for  command  and  control  inputs,  operator 
interface,  command  outputs,  and  status  reports.  An  interface  would  need  to  be  developed  for 
NILM  and  MFM  to  operate  together.  This  could  be  developed  within  the  current  NILM 
architecture.  A  typical  installation,  as  installed  on  the  Famous  Class  USCG  Cutter  (WMEC-907 
Escanaba),  is  shown  in  Figure  2. 


Figure  2:  NILM  installation  (3). 


The  NILM,  to  date,  has  been  used  in  radial  distribution  systems  and  has  not  been  explored  in  a 
ring  distribution  system  such  as  that  employed  by  DDG  IIA.  An  important  reason  for  building 
the  hardware  emulator  is  to  test  NILM  capabilities  in  a  ring  distribution  plant. 

1.3  Thesis  Objectives 

The  5000  Watt  bench  generator  shown  in  Figure  3  is  being  electrically  and  mechanically 
characterized  to  determine  feasibility  as  a  test  platform  for  NILM  fault  studies.  The  purpose  of 
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this  thesis  is  to  compare  the  characteristics  of  the  bench  generator  to  a  DDG-51  gas  turbine 
generator.  The  first  step  for  this  project  is  to  obtain  short  circuit  real  data  for  the  bench  generator 
and  then  to  match  a  simulation  to  the  real  data.  When  this  is  accomplished  it  means  that  the 
short  circuit  time  constants  as  well  as  the  transient  reactances  are  known  and  the  response  of  the 
bench  generator  can  be  predicted  for  fault  conditions.  NAVSEA  has  provided  data  necessary  to 
simulate  short  circuits  for  DDG-51  gas  turbine  generators  associated  with  3  MW  and  4.5  MW 
plants.  Since  the  reactances  and  time  constants  for  the  bench  generator  and  DDG-5 1  plants  are 
known  they  can  be  used  to  scale  and  compare  the  NILM  bench  generator  results  to  a  potential 
realization  of  NILM  implementation  on  a  DDG-51  plant.  These  results  would  need  to  be 
confirmed  in  an  actual  DDG-51  plant  and  the  goal  would  be  to  do  so  at  the  Land  Based 
Engineering  Site  (LBES)  in  Philadelphia,  PA. 


Figure  3:  Ming  Dong  5000  Watt  Bench  Generator 

The  bench  generator  is  a  4  pole,  field  wound  synchronous  generator.  The  stator  is  Y-connected 
with  the  neutral  ungrounded  and  a  rated  voltage  of  120  Vrms  phase  to  neutral.  Characterization 
of  the  generator  parameters  are  discussed  in  chapter  5.  The  bench  generator  is  being  driven  by  a 
1.5  HP  DC  motor  (120  Volts  DC,  11.2  Amps,  4800  RPM).  Leedback  loops  are  being  used  to 
control  armature  frequency  and  voltage  responses.  The  goal  is  to  relate  the  speed  and  voltage 
responses  of  the  bench  generator  to  a  DDG-5 1  plant  response. 
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Chapter  2 


The  DDG  IIA  electric  plant  was  developed  in  accordance  with  the  NGIPS  architecture  for  zonal 
electrical  distribution  systems.  The  DDG  IIA  plant  is  a  ring  distribution  system  which  means 
that  power  flowing  to  a  particular  load  may  be  coming  from  more  than  one  generator.  Having  a 
NILM  located  at  each  MFM  will  place  it  upstream  in  the  electric  plant  such  that  the  power 
flowing  to  a  particular  load  may  be  going  through  more  than  one  NILM.  Chapter  I  discussed 
that  classifying  the  bench  generator,  both  mechanically  and  electrically,  is  the  first  step  towards 
the  final  goal  of  having  a  hardware  emulator  for  the  DDG  IIA  electric  plant.  One  of  the  major 
reasons  for  the  hardware  emulator  is  to  test  the  NILM  in  a  ring  distribution  lineup  to  evaluate  if 
power  to  particular  loads  may  be  seen  from  more  than  one  NILM.  The  ability  to  identify  loads 
in  a  ring  distribution  plant  could  potentially  be  used  by  MFM  to  improve  fault  localization  times. 
This  chapter  will  provide  information  about  the  DDG  IIA  plant  and  MFM  to  show  the 
motivation  for  the  hardware  emulator  realization. 


2.0  Zonal  Electrical  Distribution  System  with  MFM  111 

Figure  4  below  shows  the  MFM  currently  employed  on  DDG-51  FLT  IIA  plants.  The  unit  is 
responsible  for  the  following  actions  (1): 

1 .  Importing  voltage,  current,  and  circuit  breaker  status. 

2.  Processing  voltage  and  current  data  to  detect  fault  conditions. 

3.  Communicate  with  the  other  MFM  throughout  the  plant. 

4.  Make  shunt  trip  decisions. 
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Figure  4:  Multi  Fuuctiou  Module  (MFM  III) 


The  MFM  III  unit  monitors  two  channels  of  data  associated  with  each  circuit  breaker.  The 
channels  correspond  to  information  on  the  upstream  and  downstream  sides  of  the  circuit  breaker. 
Each  channel  has  two  voltage  sensors  that  monitor  phase-to-phase  voltages  (Vab,  Vcb)  from 
potential  transformers  (440:110),  as  well  as  all  three  line  currents  from  current  transformers 
(6000:5).  The  MFM  units  also  have  three  ethernet  communication  ports  that  are  used  for  point 
to  point  communications  as  well  as  “ring”  communications.  All  of  the  MFM  modules  are  in 
direct  communications  with  the  MFM  that  are  adjacent  to  them.  This  is  referred  to  as  point  to 
point  communications.  These  communications  are  necessary  because  trip  logic  for  isolation 
determination  is  dependent  on  the  status  of  adjacent  circuit  breakers.  The  MFM  point  to  point 
communications  occur  every  1.0  ms.  Ring  communications  are  system  wide.  Every  5.0  ms  the 
ring  communication  network  sends  out  information  with  a  status  of  all  MFM  to  all  units.  A 
shunt  trip  signal  is  generated  based  on  each  MFM  local  voltage  and  current  measurements  as 
well  as  logic  tables  associated  with  the  other  MFM  in  the  system.  The  plant  configuration  for  the 
DDG-51  FET  II A  is  shown  in  Figure  5. 
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Figure  5:  Electric  Plant  Layout  (DDG51  FLT  IIA) 

Since  both  MFM  and  NILM  monitor  phase-to-phase  voltages  and  line  currents  there  is  an 
opportunity  to  have  both  systems  collocated  and  operating  from  the  same  set  of  inputs.  NILM 
capabilities  could  enhance  MFM  III  operation  through  optimization  of  trip  set  points  as  well  as 
adding  the  proven  benefits  of  equipment  monitoring  for  eonditions  based  maintenanee. 


2.1  High  Speed  Relay  (HSR)  Algorithm 

MFM  uses  the  HSR  algorithm  to  assess  when  fault  conditions  occur  and  to  determine  where  the 
fault  is  oeeurring  sueh  that  the  eleetric  distribution  system  only  isolates  those  affeeted  areas.  The 
algorithm  uses  the  infamous  Parks  transformation  of  the  three  phase  voltages  to  diseriminate 
faults  from  normal  ships  transients  and  uses  ehanges  in  power  flow  to  assess  where  a  fault  is 
loeated.  Figure  5  above  shows  what  the  positive  direetion  of  current  flow  whieh  is  the  same  as 
positive  power  flow.  The  Park  transformation  is  shown  below. 
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Where  the  voltage,  v,  is  transformed  from  the  rotating  referenee  frame  of  a,  b,  and  e  to  the  direct 
and  quadrature  frame  of  reference  of  d  and  q.  If  the  sum  of  Va,  Vb,  and  Vc  are  equal  to  zero  than 
Vo  will  also  equal  zero  which  is  the  assumption  used  in  the  HSR  algorithm.  The  assumption  of  a 
constant  operating  frequency  results  in  the  time  varying  functions  of  voltage  being  transformed 
into  a  constant  magnitude  and  angle  as  shown  below. 


COs(o)t  +  (p) 

(  2n 

Vfj  —  Vf,  cos  lojt  — ^  ^ 

/  2n 

Vc  =  Vc  cos  (  mt  +  ^ 


Vdq  = 


<  -cp 


The  magnitude  and  phase  of  Vdq  can  then  be  used  to  assess  whether  a  fault  condition  has 
occurred.  In  normal  plant  operations,  with  changing  load  conditions,  there  will  be  changes  in 
Vdq  but  these  will  be  a  small  percent  of  what  would  be  seen  in  a  fault  condition.  Any  changes  in 
magnitude  by  a  set  threshold  ±(Tm)  between  successive  samples  would  be  assessed  as  a  fault 
condition.  Similarly,  any  changes  in  angle,  (p,  would  indicate  changes  in  synchronization  and 
would  be  compared  to  thresholds  established  to  determine  fault  conditions. 

Once  a  fault  has  been  detected  the  HSR  algorithm  then  uses  a  running  average  of  power  to 
determine  the  direction  of  the  fault  with  respect  to  each  individual  MFM.  These  directions  are 
with  respect  to  the  direction  of  current  flow  as  shown  schematically  in  Figure  5.  Once  a  fault  has 
been  identified,  through  use  of  the  Parks  transformation  on  phase-to-phase  voltages,  power  is 
assessed  to  determine  which  breakers  need  to  be  tripped.  At  each  MFM  power  is  calculated  for 
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each  sample  and  a  running  average  is  kept.  After  a  fault  condition  is  assessed  the  average  value 
of  power,  at  each  MFM,  is  frozen  and  used  to  compare  to  successive  power  calculations  during 
the  fault.  If  power  increases  by  a  factor  of  1 .5  over  the  steady  state  average  power  at  a  particular 
MFM,  the  associated  circuit  breaker  gets  tripped.  Power  is  calculated  according  to  calculations 
shown  below. 


Power  —  ^bn^b  ^cn^c 

~  ^an^a  ^bn^b  ^cn^c  ^cnQ-a  ~  ^cnQ-b  ~  ^i>) 

~  ~  ^cn)^a  (Pbn  ~  ^cn)^b  ^cn(.^a  A) 

~  ~^ca^a  ”f  ^bc^b  ”f  ^cnO-a  ”f  ”f  A) 

=  ~^ca^a  ”f  ^bc^b 

The  power  equation  is  derived  above  to  use  the  actual  inputs  that  are  provided  from  the  MFM 
which  are  the  two  phases-to-phase  voltages  and  three  line  currents.  It  is  observed  that  this  still 
requires  a  line  to  neutral  voltage,  Vcn,  which  the  MFM  does  not  provide.  Once  again  the 
assumption  is  that  the  three  line  currents  will  sum  to  zero  in  which  case  the  line  to  neutral 
voltage  is  not  needed.  If  the  sum  of  the  currents  are  not  zero  than  the  power  calculated  is  not 
correct.  This  is  because  the  fault  can  no  longer  be  modeled  as  a  load  on  that  particular  line  and 
involves  a  connection  to  another  point  in  the  system  (2).  In  this  case  more  than  one  MFM  would 
activate  a  shunt  trip  in  order  to  fully  isolate  the  fault. 
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Chapter  3 


The  previous  chapters  have  discussed  the  motivation  for  building  an  emulator  for  the  DDG  IIA 
plant.  The  first  step  towards  building  an  emulator  is  to  characterize  the  bench  generator  that  will 
be  used.  To  accomplish  this,  laboratory  experiments  will  be  conducted  to  derive  electrical 
constants.  Once  identified,  the  electrical  constants  can  be  used  to  simulate  the  bench  generator. 
This  chapter  will  discuss  the  background  of  the  electrical  generator  that  is  being  evaluated  for 
use  in  the  hardware  emulator.  Also  discussed  is  the  simulation  that  will  be  used  for  the  bench 
generator  and  the  DDG  generator. 

3.0  Synchronous  Machines 
3.1  Basic  Machine  Layout 

Shown  in  Figure  6  is  a  three  phase,  four  pole,  Y-connected  synchronous  generator.  The  bench 
generator  is  a  salient  poled  generator.  The  armature  winding  is  on  the  stator  and  the  field 
winding  is  on  the  rotor.  The  field  winding  is  excited  by  a  DC  voltage  source  which  is  applied 


Figure  6:  Four  pole,  Y-Couuected  Syuchrouous  Geuerator  (5) 


through  a  set  of  stationary  brushes  which  are  in  contact  with  rotating  slip  rings.  Since  the 
generator  has  two  sets  of  poles,  every  time  the  rotor  rotates  one  full  revolution  the  generator 
electrical  frequency  rotates  two  full  revolutions.  This  results  in  the  generator  having  an  electrical 
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frequency  of  60  Hz  when  the  mechanical  frequency,  n,  of  the  generator  is  1800  RPM  as  shown 
in  the  equation  below. 


3.2  Direct  and  Quadrature  Transformation 


The  synchronous  generator  can  be  represented  by  6  windings.  The  windings  are  the  three 
armature  windings,  the  field  winding,  and  two  damper  windings  (one  on  the  direct  axis  of  the 
field  winding  and  the  other  on  the  quadrature  axis).  The  equations  representing  the  synchronous 
machine  are  as  shown  below  (6). 


—Ph 

hhK' 

Lph 
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-Lr  - 
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^ab  ^b 

^bc 

Lr  — 

Lfkd 

Lkd 

0 

-^ac  ^bc 

Lei 

0 

0 

Lkq 

From  these  equations  it  can  be  seen  that  the  direct  axis  rotor  flux  is  not  dependent  on  the 
quadrature  damping  currents.  The  same  is  true  that  the  quadrature  axis  flux  does  not  depend  on 
the  field  current  or  the  direct  axis  damping  current. 
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M  = 


The  stator  to  rotor  mutual  inductance  (  M  )  relates  the  phase  fluxes  to  the  rotor  currents  and  the 

rotor  fluxes  to  the  phase  currents.  The  stator  to  rotor  mutual  inductance  is  dependent  on  the 
angle  of  the  rotor  to  the  line  of  flux  aligned  with  the  phase  A  current  and  as  such  will  have  a  time 
dependency  as  the  rotor  rotates. 


ifl  =  iflo  +^2  cos  2(0) 
h  =  iflo  +  ^2  cos  2  (^0  -  y) 

=  ^ao  +  h  cos  2(0+^^ 

^ab  ~  ^abO  ^2  COS  2  ^0  —  — ^ 

^bc  ^abO  ^2  COS  2(0) 

^ac  ~  ^abO  ^2  COS  2  ^0  +  — ^ 

For  the  bench  generator  there  is  saliency  which  explains  the  form  of  the  equations  above  for  La, 
Lb,  Lc,  Lab,  Lbc,  and  Lac.  L2  is  the  same  in  each  of  the  expressions  and  0  is  the  position  of  the 
rotor  relative  to  phase  A.  The  Parks  transformation  will  be  used  to  transform  to  a  coordinate 
system  in  which  the  rotor  is  stationary.  For  the  proof  shown  below  T_  is  the  Parks  transformation 

and  is  the  inverse  Parks  transformation. 

^ph  =  Lphiph  +  KLr 

Now  apply  the  Parks  transformation.  In  the  equation  below  we  set  Lph  —  r  ^Ldq- 
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L^ph-LyjhL  ^Ldq  +ZKLr 
Since  T_  Apf^  =  we  ean  rewrite  the  equation  as  shown  below. 


Zdq  ZdqLdq  ”1”  Zc  Lr 

Where, 

Ldq  —  L  ^ 

Lc  =  TM 


Using  the  same  methodology  one  ean  obtain  the  following  equations  for  the  rotor  fluxes. 

ZR=K^L~^Ldq+^LR 

_3  ^ 

“  2ZcLdq  Lr  Lr 

Where, 

Putting  it  altogether  in  one  equation  we  have. 


Ldq  and  Lc  have  the  eonvenient  expressions  shown  below  (6). 
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This  is  then  separated  into  three  current  flux  relationships  for  the  direct  axis,  indirect  axis,  and 
zero  sequence  components.  This  can  be  accomplished  because  these  axes  are  orthogonal  in 
space  (6).  The  zero  sequence  component  is  neglected  and  the  direct  and  quadrature  relationships 
are  shown  below.  The  zero  sequence  component  can  be  neglected  because  any  coupling  to  the 
direct  axis  is  through  higher  harmonics  that  do  not  couple  well  to  the  armature  (6). 
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3 

Id' 

^kd 

= 
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Equation  1:  Direct  and  Quadrature  Flux  Current  Relationship 


Now  the  above  results  can  be  used  to  solve  for  the  direct  and  quadrature  voltages. 


Therefore, 


Where, 
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d.  d. 
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Which  results  in  the  final  equation  for  direct  and  quadrature  voltages  shown  below. 


dt 
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dAa 
^  dt 


(jOAn 


It  should  be  noted  that  the  voltages  listed  above  do  not  inelude  the  effects  of  the  armature 
resistance.  To  obtain  the  terminal  voltage  the  armature  resistance  would  need  to  be  accounted 
for  which  will  be  shown  in  the  next  section. 


3.3  Equal  Mutuals  Model  Description 

Equation  1 ,  derived  in  the  previous  section,  is  an  equation  in  which  the  inductance  matrix  has  six 
variables.  These  six  variables  are  unknown  for  the  bench  generator.  The  equal  mutuals  model  is 
a  model  that  sets  the  off  diagonal  terms  equal  to  each  other  which  allows  a  circuit  model  to  be 
determined(6).  From  this  simplification  it  is  possible  to  determine  the  characteristics  of  the 
bench  generator  that  will  allow  a  simulation  to  be  developed. 

Equation  1  in  the  previous  section  will  now  be  written  in  per  unit  form  as  shown  below.  A  direct 
comparison  of  the  electrical  parameters  of  the  DDG  and  bench  generator  is  difficult  since  the 
machine  ratings  are  different.  Per  unit  values  will  be  used  to  allow  a  direct  comparison  of  the 
electrical  constants  for  the  bench  and  DDG  generators. 


^d 

■  Xa 

^akd 

^ad 

id' 

4’kd 

^akd 

^kd 

^fkd 

ikd 

.  ^ad 

^fkd 

[if\ 

%  ' 

Xq  Xakq'\ 

iq 

^kq- 

.  ^akq  ^kq  \ 

ikq. 

The  bases  have  been  chosen  such  that  the  inductance  matrices  are  symmetric  (6).  The  “equal 
mutuals”  model  will  now  be  used  to  set  the  field  and  damper  base  impedances  such  that  all  three 
mutual  inductances  are  equal  (6). 


^akd  ^fkd  ^ad 
^aq  ~  ^akq 
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Where  Xakq  is  now  to  be  referred  to  as  Xaq  so  that  the  quadrature  model  is  consistent  with  the 
direct  axis  model. 


Figure  7:  Direct  Axis  Model 


ra  Xal 


Figure  8:  Quadrature  Axis  Model 
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Figure  7  and  Figure  8  above  show  the  physieal  model  assoeiated  with  Equation  2  with  the 
leakage  reaetances  defined  as  follows.  (6) 


^al  ~  ^ad 

^kdl  ~  ^kd  ~  ^ad 
^fl  ~  ~  ^ad 

With  this  model  the  same  problem  as  before  exists  which  is  how  to  obtain  the  relevant 
parameters  to  be  able  to  evaluate  the  synchronous  generator.  This  is  discussed  in  chapter  five 
but  it  turns  out  that  results  from  short  circuit  testing  and  modeling  can  be  used  to  obtain  the 
generator  parameters. 
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Chapter  4 


This  chapter  will  discuss  all  of  the  laboratory  testing  and  results  from  the  bench  generator. 

4.0  Preliminary  Data  Collection 

To  be  able  to  compare  the  bench  generators  to  any  Navy  electrical  generator  it  is  first  necessary 
to  classify  the  characteristics  of  the  bench  generator.  The  preliminary  tests  accomplished  in 
order  to  characterize  the  bench  generator  were  the  following: 

1 .  DC  test:  Armature  Resistance  (Ri).  (7) 

2.  Open  Circuit  Test:  The  Magnetization  Characteristic  and  Field  Resistance.  (7) 

3.  Short  Circuit  Test:  Synchronous  D-Axis  Reactance  (Xd).  (7) 

The  specifications  of  the  bench  generator  are  shown  below  in  Table  1.  This  generator  is  a  salient 
pole  synchronous  wound  field  machine. 


Power  (W) 

5000 

Voltage  (phase  to  phase),  (Vrms) 

207.85 

Current  (phase  to  neutral),  (Arms) 

13.89 

Frequency  (Hz) 

60 

Poles 

4 

Power  factor 

0.8 

Table  1:  Bench  Generator  Specifications 


4.1  DC  test 

The  bench  generator  is  a  Y-connected  generator,  with  the  neutral  ungrounded,  as  pictured  below 
in  Figure  9.  As  a  result  the  armature  resistance  (Ri)  is  obtained  simply  by  averaging  the  three 
armature  resistances  (Ra,  Rb,  and  Rc)  and  multiplying  by  temperature  and  AC  skin  effect 
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correction  factors  (Ktemp  and  Kac).  A  typical  value  for  Ktemp  multiplied  by  Kac  is  1.5  (7).  Lab 
measurements  yielded  an  armature  resistanee  equal  to  0.48  Q.. 

P,  _  ^ca 

Rdc-  g 

Now  to  account  for  temperature  and  the  AC  skin  effect  Rdc  is  multiplied  by  1.5. 

i?l  =  KjEMpRacRdc 


Figure  9:  Bench  Generator  Y-configuration(7) 


4.2  Open-Circuit  Test:  Magnetization  Characteristic  and  Field  Resistance 

The  Open  Cireuit  Charaeteristies  (OCC)  are  obtained  with  the  generator  at  synehronous  speed  by 
measuring  the  open  eircuited  armature  phase  to  phase  voltage  (Etpjp)  while  the  field  current  (Ip)  is 
varied  such  that  the  terminal  voltage  goes  from  zero  to  rated  voltage.  While  Ip  is  being  adjusted 
the  field  voltage  (Vp)  is  also  being  measured  and  therefore  the  field  resistance  (Rp  =  Vp/Ip)  can 
be  averaged  throughout  the  range  of  operation.  The  average  Rp  was  equal  to  15.8  O. 

Using  the  equation  shown  below  and  the  open  cireuit  eharaeteristics  lab  results  we  ean  evaluate 
Kag,  Ax,  and  Bx. 
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Ef  „  p 

Eag 

Where  Ax  is  the  change  in  field  current  which  accounts  for  saturation,  Bx  is  an  exponential 
constant  that  determinations  the  exponential  shape  of  the  field  current,  and  Kag  slope  of  the  air 
gap  line  (linear  portion  of  curve). 

The  results  are  shown  below  in  Figure  10.  The  results  were  as  expected  and  we  see  the  stator 
voltage  beginning  to  saturate  as  If  is  increased. 
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Figure  10:  Bench  Generator  OCC 
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4.3  Short  Circuit  Characteristics 


The  OCC  results  from  the  previous  section  along  with  the  short  circuit  characteristics  (SCC), 
which  will  be  determined  in  this  section,  will  be  used  to  determine  the  direct  axis  reactance. 
Using  the  setup  shown  in  Figure  1 1  the  synchronous  generator  was  operated  at  rated  frequency 
while  the  stator  terminals  were  short  circuited.  If  was  slowly  adjusted  from  zero  amps  up  to  741 
mA  while  the  Ia  was  measured.  U  and  Ip  shown  in  Figure  1 1  are  meters  used  to  measure  the 
field  and  stator  currents. 


Figure  11:  SCC  and  OCC  test  setup  (7) 

The  results,  called  the  short  circuit  characteristics  (SCC),  are  linear  since  the  air-gap  flux  is  low 
and  the  stator  is  magnetically  unsaturated  at  low  levels  of  Ip.  The  results  are  shown  in  Figure  12 
and  Figure  13. 
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Figure  12:  Bench  Generator  Iscc 
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Figure  13:  Bench  Generator  OCC  and  SCC 

The  next  step  in  the  process  is  taking  the  OCC  and  SCC  and  obtaining  the  d-axis  reactance. 
Under  short  circuit  conditions,  the  stator  current  lags  the  armature  voltage  by  nearly  90*^  which 
results  in  Ia  being  approximately  100%  U.  Therefore,  Xd  is  obtained  by  the  following  equation: 
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^  _  Efjlp  =  Ipo)fromOCC 
^  “  IaOF^lFo)fromSCC 

=  7^(4  =  13.89i4)  from  SCC  =  741mA 

One  can  see  that  when  If  is  relatively  small  that  the  armature  voltage  will  be  approximately 
linear.  The  air  gap  reactance  Xdag  is  defined  for  small  levels  of  If  and  is  relatively  constant. 
However,  as  the  stator  becomes  saturated  Xd  will  begin  to  decrease.  The  variable  Xd  is  defined 
as  that  value  of  reactance  that  corresponds  to  rated  armature  current  in  Figure  13  above  and  as 
seen  in  Figure  14  below  corresponds  to  8.70. 


9.5 

9 

8.5 

8 

X°  7.5 
7 

6.5 

6 

5  5 

■  0  0.5  1  1.5  2  2.5 

'f 

Figure  14:  Bench  Generator  D-Axis  Reactance 
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4.4  Short  Circuit  Testing 


The  short  circuit  test  requires  the  bench  generator  to  be  operating  at  rated  frequency  and  voltage 
open  circuited.  A  switch  is  then  used  to  short  the  three  armature  terminals  simultaneously.  This 
is  not  to  be  confused  with  the  SCC  testing  that  was  accomplished  in  the  previous  section.  The 
parameters  that  are  obtained  from  this  test  are  the  transient  and  sub-transient  short  circuit  d-axis 
time  constants  (T^  and  T^)  as  well  as  the  transient  and  sub-transient  d-axis  reactances 
and  X^).  The  method  employed  to  determine  the  constants  was  to  use  the  equation  shown 
below  which  is  called  the  “classical  short  circuit”  equation(8). 

^\Eam\—+(^-—]e  +  'Ad]sin(mot  +  5) 

Equation  3:  Classical  Short  Circuit  Equation 

Where  Ea(0)  is  the  line  to  neutral  rated  voltage  which  is  120  V  for  the  generator  and  the 
sinusoidal  term  can  correspond  to  any  of  the  three  phases.  The  importance  of  this  equation  is  the 
exponential  envelope.  There  are  four  degrees  of  freedom  in  this  equation  (A^,  A^',  T^,  and  T^) 
so  an  assumption  is  necessary  to  obtain  the  machine  data.  The  assumption  that  was  used  was 
that  the  dominant  reaction,  in  fact  maybe  the  only  reaction,  would  be  as  a  result  of  the  transient 
parameters.  At  the  start  of  testing  for  the  5000  Watt  generator  there  was  no  initial  premise  that 
any  sub-transient  reaction  would  be  seen  after  looking  at  the  short  circuit  data. 
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Figure  15:  Short  circuit  euvelope  (classical  equatiou)  (8) 

Figure  15  above  shows  three  regions  of  short  eireuit  eurrent.  The  first  region  is  the  sub-transient 
region  and  the  envelope  is  defined  by  X'J^  and  .  The  seeond  region  is  the  transient  region  and 
is  defined  by  and  T^.  Finally,  steady  state  is  purely  a  funetion  of  the  d-axis  reaetanee  (Xa). 
By  matching  the  short  circuit  currents  obtained  with  laboratory  testing  it  is  possible  to  determine 
all  short  circuit  reactances  and  time  constants. 


4.4.1  Test  Setup 

The  DC  motor  shown  in  Figure  16  was  used  for  the  initial  testing  of  the  generator.  The  motor  is 
only  1.5  HP  which  would  equal  1119  Watts  of  electrical  power.  The  DC  motor  would  not  be 
appropriate  to  be  the  prime  mover  for  the  bench  generator  if  we  were  operating  under  full  load 
conditions  but  for  purposes  of  obtaining  constants  for  the  generator  it  was  an  adequate  motor. 
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Figure  16:  Bench  Generator  Prime  Mover  (1.5  HP) 

The  motor  could  easily  maintain  the  generator  at  60  Hz  but  slowed  down  excessively  during  any 
loading  conditions.  This  was  especially  prevalent  during  the  short  circuit  testing.  A  tachometer 
was  used  to  capture  the  speed  profile  of  any  transients  which  was  necessary  to  be  incorporated 
into  simulations  later  on.  The  switches  shown  in  Figure  17  were  used  to  insert  loading 
conditions  as  well  as  short  circuits  in  a  controlled  manner. 


Figure  17:  Switching  Apparatus  for  short  circuit  implementation 
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TTL  logic  was  used  to  have  an  adjustable  trigger  switch.  A  switch  on  duration  of  52  ms  was 
used  when  the  generator  was  being  short  circuited.  The  generator  was  always  open  circuited  in 
the  initial  condition  and  loads  of  progressively  smaller  resistances  were  applied  to  the  generator. 
In  much  the  same  manner  as  we  approached  the  fault  condition  the  on  times  of  the  switch  were 
raised  until  the  generator  reached  52  ms. 


4.4.2  Test  Results 

Figure  18  below  shows  typical  lab  results  obtained  for  the  bench  generator  short  circuit  tests. 
These  results  show  the  tachometer  results  for  the  1 .5  HP  motor.  The  tachometer  is  a  Hall  effect 
sensor  that  was  placed  in  close  proximity  to  the  gears  on  the  uncoupled  side  of  the  motor.  The 
output  of  the  tachometer  was  used  to  capture  the  speed  profde  of  the  generator  during  the  short 
circuit  transient. 
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Figure  18:  Short  Circuit  Results  from  Oscilloscope 
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The  TTL  voltage  input  to  each  of  the  three  switches  is  also  shown.  In  this  example  it  was  noted 
that  there  was  a  6.6  ms  delay  from  the  time  the  voltage  was  applied  to  the  switch  until  the  switch 
was  shut  as  indicated  by  the  delay  in  short  circuit  line  current.  The  same  delay  was  apparent  in 
the  opening  of  the  switch.  The  delays  were  consistent  for  each  of  the  three  switches  and  did  not 
affect  the  short  circuit  results. 


Figure  19  below  shows  the  field  current  response  during  the  short  circuit.  In  the  steady  state 
prior  to  the  short  circuit  the  field  current  is  constant  (approximately  two  amps)  and, 
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Figure  19:  Field  Current  during  short  circuit 


correspondingly,  the  flux  linkages  in  the  field  winding  are  constant.  At  this  point  the  stator 
currents  are  zero  and  therefore  do  not  contribute  to  the  field  flux  linkages.  Once  the  stator 
currents  start  in  the  short  circuit  condition  they  lag  the  internal  voltages  by  90^  which  causes  a 
demagnetizing  armature  reaction  flux  in  the  field  winding.  This  armature  reaction  opposes  the 
field  flux.  From  the  conservation  of  flux  linkages,  and  assuming  the  field  resistance  equal  to 
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zero  we  get  the  result  that  the  field  fiux  will  remain  constant.  The  field  resistance  is  not  zero,  but 
is  small  in  value,  therefore  for  a  small  time  interval  such  as  seen  above  we  can  assume  that  the 
field  flux  will  be  approximately  constant  (8).  As  a  result  of  the  demagnetizing  component  of  the 
stator  currents  the  field  current  must  increase.  Figure  20  shows  the  speed  profile  of  the  bench 
generator.  The  first  speed  profile  shown  is  as  a  result  of  the  tachometer  connected  to  the  1.5  HP 
DC  motor.  The  second  profile  is  as  a  result  of  the  simulation  which  will  be  discussed  in  the  next 
chapter.  The  results  match  each  other  with  the  exception  of  the  choppiness  of  the  tachometer 
waveform.  For  the  simulation  the  speed  profile,  based  on  the  inertia  constant  was  used.  The 
parameter  necessary  for  converting  the  tachometer  to  frequency  is  the  gear  ratio  of  the  generator 
to  the  motor  which  is  (12/32).  When  the  DC  motor  is  operating  at  4800  RPM  the  generator  is 
turning  at  1800  RPM  which  equates  to  an  electrical  frequency  of  60  Hz. 


Figure  20:  Speed  profile  for  Short  Circuit  (Tachometer  aud  luertia  Coustaut) 
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4.4.3  Analysis  of  Data 


As  stated  previously,  the  purpose  of  the  short  circuit  testing  was  to  obtain  the  transient  and  sub¬ 
transient  reactances  and  X'J)  as  well  as  the  transient  and  sub-transient  short  circuit  time 

constants  and  t‘J)  for  the  bench  generator.  Figure  21  shows  how  was  evaluated. 


Figure  21:  Evaluation  of  Xd' 

There  is  a  typical  range  of  values  (9)  associated  with  direct  axis  transient  reactances  (0.17-0.37 
per  unit)  which  was  the  starting  point  for  the  determination  of  the  short  circuit  time  constants. 
The  reactances  must  be  matched  to  Figure  21  above,  matching  the  classical  short  circuit 
equation,  as  well  the  bench  generator  “Equal  Mutuals”  model  which  will  be  discussed  in  the  next 
chapter.  The  method  employed  was  to  assume  that  there  was  no  sub-transient  reaction  and  that 
the  short  circuit  response  was  entirely  as  a  result  of  transient  reactance.  The  transient  envelope 
using  Equation  3  is  shown  in  Eigure  21  above.  The  short  circuit  current  corresponding  to  60  ms 
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was  matched  to  the  transient  envelope  as  this  was  the  last  actual  datum  point  available.  It  is  then 
apparent  that  the  initial  assumption  of  no  sub-transient  component  is  not  a  valid  one.  There  is  a 
sub-transient  component  for  the  generator  response  since  the  short  circuit  current  is  not  bound  by 
the  transient  envelope.  Therefore,  what  is  shown  above  is  the  transient  envelope  that  was 
matched  to  the  “Equal  Mutuals”  model.  This  corresponded  to  a  short  circuit  transient  time 
constant  of  45  ms  and  a  transient  reactance  of  0.197  Q.  per  unit. 


Figure  22:  Xd"  Determination 

It  is  now  possible  to  obtain  the  sub-transient  reactance  and  short  circuit  time  constant.  In  a 
similar  fashion  to  choosing  the  transient  reactance,  X'J  was  chosen  to  satisfy  the  “classical”  short 
circuit  equation  as  shown  in  Figure  22  above  as  well  as  producing  a  good  result  in  the  chapter 
five  simulations.  T^'  can  then  be  chosen  to  match  the  short  circuit  response  to  the  sub-transient 
envelope  shown  above.  There  is  not  a  perfect  match  of  the  sub-transient  envelope  to  the  short 
circuit  current  shown  above.  The  solution  of  T^'  equal  to  15  ms  and  equal  to  0.094  Q.  per 
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unit  was  chosen  based  on  minimizing  errors  with  both  the  elassieal  equation  as  well  as  the 
simulation  discussed  in  the  next  ehapter. 
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Chapter  5 


This  chapter  takes  the  laboratory  results  from  ehapter  4  and  uses  the  equal  mutual  model  to 
simulate  the  DDG  and  beneh  generators.  The  simulation  models  allow  direet  eomparison  of  the 
two  generators.  This  ehapter  will  summarize  those  results  and  determine  feasibility  for  hardware 
emulation. 

5.0  Generator  Modeling 

Now  that  the  physieal  parameters  of  the  beneh  generator  are  known  the  next  step  is  to  use  those 
parameters  to  simulate  the  generators  short  eireuit  response.  Onee  the  simulation  matehes  the 
actual  response  for  the  bench  generator  it  is  then  possible  to  use  the  DDG-5 1  plant  parameters  to 
simulate  its  short  eireuit  response.  At  this  point  we  ean  eompare  the  differences  between  the  two 
generator  responses  and  present  a  eonclusion  on  how  the  beneh  generator  results  ean  be  used  to 
assess  how  the  DDG-5 1  plant  would  respond.  Validation  of  the  model  used  will  be 
aeeomplished  by  looking  at  aetual  short  eireuit  data  (10).  The  simulation  results  are  assumed 
validated  when  the  simulation  results  mateh  the  known  response. 

5.1  Equal  Mutuals  Modeling 

The  list  of  parameters  that  were  unknown  prior  to  the  short  eireuit  analysis  are  shown  below. 
These  are  the  values  that  are  neeessary  to  be  identified  in  order  to  aeeomplish  a  short  eireuit 
simulation  of  the  synehronous  generator. 

1 .  Inertia  Constant  H  -  This  eonstant  relates  the  aetual  moment  of  inertia  of  the  meehanieal 
rotating  eomponents  to  a  base  of  eleetrieal  volt  amps.  The  speed  profile  slowdown  is 
known  for  the  generator  as  a  result  of  the  taehometer.  The  inertia  eonstant  is  used  to 
mateh  the  model  speed  ehange  to  that  seen  in  the  laboratory. 

_  Rotational  kinetic  energy  at  rated  speed  Jcoq 
Base  Power 
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2.  Direct  Axis  Reactances  and  A^')  -  These  values  have  a  direct  effect  on  both  the 
Equal  Mutuals  model  and  with  the  short  circuit  “classical  equation”  analysis  in  the 
previous  chapter.  These  parameters  are  obtained  by  matching  the  simulation  model  and 
the  classical  equation.  The  short  circuit  time  constants  can  then  be  determined  from  the 
classical  equation  with  the  assumption  that  the  dominant  reaction  for  the  bench  generator 
would  be  from  the  transient  reaction.  The  remainder  of  the  short  circuit  reaction  would 
be  attributed  to  the  sub-transient  reaction. 

3.  Short  Circuit  Direct  Axis  Time  Constants  (T^  and  T^)  -  These  constants  determine  the 
amplitude  of  the  short  circuit  current  response  in  the  generator.  Based  on  the  analysis  in 
Chapter  four  we  have  an  initial  range  of  values.  The  short  circuit  time  constants  do  not 
have  a  direct  impact  on  the  equal  mutuals  model  although  they  can  be  used  to  assess 
initially  what  other  parameters  could  be. 

4.  Open  Circuit  Direct  Axis  Time  Constants  (T^g  and  -  When  the  armature  is  open 
circuited  the  field  circuit  is  no  longer  affected  by  the  armature.  The  change  of  field 
current  is  dependent  only  on  the  field  winding  self-inductance  and  resistance.  Therefore, 
the  open  circuit  armature  voltage  is  proportional  to  field  current  and  has  the  same  time 
constant  as  the  field  circuit.  These  values  are  determined  from  the  equal  mutuals 
simulation. 

5.  Quadrature  Axis  Reactances  (A^andAg)  -  The  definitions  for  the  quadrature  axis 
reactances  are  the  same  as  for  the  direct  axis  reactances  except  that  the  rotor  is  positioned 
such  that  the  quadrature  is  in  line  with  the  crest  of  the  magneto  motive  force  wave.  For 
the  simplified  two  pole  case  shown  in  Figure  23  below  it  is  apparent  that  the  inductance 
for  the  quadrature  axis  would  be  less  than  the  direct  axis  inductance.  This  is  not, 
necessarily,  the  case  for  the  sub-transient  reactances  and  in  fact  the  q-axis  sub-transient 
reactance  ends  up  being  greater  than  the  d-axis  sub-transient  reactance  for  the  bench 
generator. 
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d  axis 


Figure  23:  Direct  and  Quadrature  Axes  (6) 


6.  Open  Circuit  Quadrature  Axis  Sub-Transient  Time  Constant  (T^'o)  -  In  a  salient  pole 
machine  the  transient  time  constant  for  the  q-axis  can  be  neglected  (11)  but  must  be 
included  in  the  model.  Although  the  generator  does  not  have  explicit  amortisseur 
windings  there  are  current  paths  that  produce  the  same  effect.  This  parameter  is 
determined  in  the  model. 

7.  Armature  Leakage  Reactance  (Xi)  -  This  is  the  reactance  which  results  from  the 
difference  between  the  total  flux  produced  by  the  armature  current  acting  alone  and  the 
space  fundamental  of  the  air  gap  flux.  There  are  three  components  that  make  up  Xi 
which  are  slot  leakage,  end  winding,  and  differential  leakage  flux.  This  parameter  is 
determined  from  the  model. 
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5.1.1  Bench  Generator  Model 


Referring  to  Figure  7  the  equations  for  the  direet  axis  reaetanee  as  well  as  the  direct  axis 
transient  and  sub-transient  reactances  can  be  written  as  shown  below. 


~  ^al  ^ad 


The  direct  axis  reactance  is  dependent  on  the  magnetizing  reactance  as  seen  in  Figure  7. 


X'd  =  Xai  +  Xad  Xfi 


The  transient  reactance  includes  the  effects  of  the  field  winding  reactance. 


Xd  ~  Xal  T  Xg^d  Xj^iWXj^di 


Finally,  the  sub-transient  reactance  is  determined  by  also  including  the  direct  axis  damping 
reactance.  These  are  determined  by  circuit  analysis  of  the  circuit  model  shown  in  Figure  7. 

These  equations  show  that  the  transient  reaction  is  defined  by  the  field  and  magnetizing 
reactances  in  parallel  whereas  the  sub-transient  reaction  includes  the  effects  of  the  damping 
reactance.  From  these  equations  we  can  solve  for  the  Xf  and  Xkd  (see  Appendix  3  and  the 
MATLAB  function  MI.M  for  equations).  The  equations  listed  below  can  be  used  to  solve  for  the 
field  and  damping  resistances.  However,  it  is  important  to  note  that  these  equations  depend  on 
the  time  constants  being  spread  widely  apart  (6).  It  should  also  be  noted  that  the  bench 
generators  field  resistance  was  directly  measured  as  opposed  to  using  a  calculated  value. 


Once  again  by  looking  at  the  circuit  model  shown  in  Figure  7  and  knowing  that  the  transient  time 
constant  is  a  function  of  the  field  winding  and  the  sub-transient  time  constant  includes  the  effects 
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of  the  damper  winding  it  is  determined  that  the  sub-transient  open  circuit  time  constant  is  the 
following; 


T'll  _ 

do  — 


^kdl  ^fl  ll^ad 

(xioTf 


Initial  estimates  for  T'do  and  T^'g  can  be  obtained  by  using  the  relationships  shown  below  which 
are  derived  from  Equation  5  with  the  assumption  that  is  much  greater  than  and  that  7’do  is 
much  greater  than  T^'g.  This  turns  out  to  be  an  adequate  assumption  to  obtain  an  initial  estimate 
for  both  andr^o-  (6) 


7’do  = 


T' f  f  _ 

do  — 


T'a 


The  equations  showed  to  this  point  allowed  for  the  calculation  of  the  parameters  that  define  the 
bench  generator  which  are  shown  in  Table  2.  This  was  not  finalized  on  the  first  iteration.  This 
required  several  iterations  which  required  looking  at  the  simulated  short  circuit  response  and 
comparing  it  to  the  actual  response  as  well  as  the  calculated  response  using  the  classical  short 
circuit  equation.  Once  the  simulated  short  circuit  response  matched  the  actual  response  those 
were  the  values  that  were  assessed  as  correct  and  are  as  listed  below. 


Phase  (Watts) 

5000 

Xkdl(f^)  p.U. 

.0628 

V base  (V)iine-line  RMS 

207.8 

Xq(a)  p.U. 

.770 

Ibase  (A)r]mS 

13.9 

Xd(Q)  p.U. 

1.007 

Pf 

0.8 

Xkq(f2)  p.U. 

.837 

Pole  Pairs 

2 

Xq  (fi)  p.U. 

.150 

H  (sec) 

.190 

Xd  (Q)  p.U. 

.197 

Ta  (Q)  p.U, 

.056 

Xd  ”(fi)  p.U. 

.094 

X|  (Q)  p.U. 

.0500 

xi<a)  p.U. 

1.13 

Xa,  (Q)  p.U. 

.720  ^  Xkd(f^)  p.U. 

1.020 
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Tkq  (Q)  p.U. 

.0317 

tdo'  (sec) 

.409 

Xkql  (ft)  p.U, 

.116 

tqo”  (sec) 

.070 

Xad  (ft)  p.U. 

.957 

tdo”  (sec) 

.030 

rf(ft)  p.U. 

.0126 

td  (sec) 

.045 

Xfl(ft)  p.U. 

.173 

td  (sec) 

.015 

rkd(ft)  p.U, 

.0185 

Table  2:  Bench  Generator  Parameters 


The  values  shown  in  yellow  in  Table  2  are  the  parameters  that  were  adjusted  in  the  simulation  to 
match  the  simulated  short  circuit  response  to  the  actual  short  circuit  response.  These  parameters 
had  a  direct  effect  on  other  calculated  values  shown  in  the  table.  The  inertia  constant  was  also 
changed  but  this  was  just  to  match  the  simulated  frequency  to  the  actual  short  circuit  frequency. 
The  inertia  constant  had  no  effect  on  the  other  calculated  values  in  the  table. 

The  model  state  equations  shown  below  are  derived  with  the  assumption  that  the  terminals  are 
constrained  by  voltage  (6).  These  are  the  equations  used  to  model  the  bench  generator  and 
DDG-51  plants. 


drpd 

dt 


=  OJoVd  +  a)xpq  -  0)oVd 


This  shows  that  the  time  rate  of  change  of  direct  axis  flux  is  dependent  on  direct  axis  voltage. 
The  base  frequency  ((Uq)  for  this  machine  is  377  rad/sec  and  does  not  change  in  the  simulation 
model.  However,  the  electrical  frequency  of  the  machine  (m)  will  slow  down  during  the  short 
circuit  scenario.  The  direct  axis  voltage  when  the  generator  is  short  circuited  is  zero  which 
simplifies  the  differential  equation. 

=  (DoVq  -  (DXpd  -  (^oralq 
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The  time  rate  of  ehange  of  quadrature  axis  flux  is  eomparable  to  the  direct  axis  flux  with 
quadrature  and  direct  subscripts  switched.  Similarly,  in  a  short  circuit  scenario  the  quadrature 
voltage  is  zero. 


dt 


—  ^O^kd^kd 


^^kq 

dt 


^O^kq^kq 


The  time  rate  of  change  of  direct  and  quadrature  damping  fluxes  are  dependent  simply  on  their 
respective  resistances  and  transient  currents.  The  resistances  are  the  same  as  those  shown  in  the 
circuit  models  (Figure  7  and  Figure  8). 


dijjf 

dt 


= 


The  time  rate  of  change  of  the  field  flux  depends  on  the  field  voltage  and  field  current.  In  the 
model  used  the  field  voltage  was  set  to  a  constant  voltage  (32  volts)  prior  to  the  short  circuit  test. 
This  field  voltage  when  the  armature  is  open  circuited  corresponded  to  rated  armature  voltage  of 
208  volts  RMS  phase  to  phase.  The  appolied  field  voltage  remains  constant  during  the  transient. 
The  actual  field  current  during  the  transient  fault  was  inserted  into  the  simulation  as  shown  in  the 
MATLAB  function  SF2.M  in  Appendix  3. 


do) 

dt 


(jOn 


Te  ~  ^d^q  ^q^d 

Equation  4:  State  Equations  for  Equal  Mutuals  Model 

The  slowdown  of  the  bench  generator  was  captured  by  using  a  tachometer.  The  tachometer 
results  were  then  used  to  identify  the  inertia  constant  for  the  bench  generator  (H)  shown  in  the 
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equation  above.  Sinee  the  generator  was  initially  open  eireuited  the  meehanieal  torque  Tm  is 
zero.  The  inertia  eonstant  was  ehanged  until  the  simulation  speed  profile  matehed  the 
taehometer  speed  profile.  The  simulation  results  were  used  instead  of  the  taehometer  speed 
profile  sinee  they  matehed. 

It  should  be  noted  that  for  short  eireuit  eonditions  the  direet  and  quadrature  voltages  (vd  and  Vq) 
are  zero  whieh  will  simplify  the  equations  above.  The  armature  resistanee  (ra),  field  resistanee 
(rf),  field  eurrent  (if),  and  speed  profile  (ro)  were  all  direetly  measured  in  the  laboratory  for  the 
beneh  generator.  All  were  direetly  input  into  the  model  with  the  exeeption  of  the  speed  profile. 
The  speed  profile  was  matehed  with  the  inertia  eonstant,  as  diseussed  previously,  and  then 
entered  into  the  model.  This  was  aeeomplished  to  determine  the  inertia  eonstant  for  the 
generator.  For  the  short  eireuit  seenario  the  torque  angle  (5)  is  initially  zero  when  open  eireuited 
and  remains  zero  during  the  short  eireuit  transient.  As  a  result  this  equation  is  not  used  in  the 
short  eireuit  analysis.  It  is  shown  beeause  it  was  used  in  the  validation  of  the  model  as  shown  in 
Appendix  6. 

The  MATLAB  simulation  model  (Appendix  3)  results  are  shown  in  Figure  24  below  for  three 
phase  short  eireuit  eurrents.  The  simulation  model  (Appendix  3)  has  three  funetions  that  are 
ealled  whieh  are  MI.M,  MIQ.M,  and  SF2.M.  The  funetion  MI.M  ealeulates  the  direet  axis 
parameters  while  MIQ.M  ealeulates  the  quadrature  axis  parameters.  The  quadrature  axis 
parameters  are  ealeulated  in  the  same  way  as  the  direet  axis  parameters  aeeording  to  the  eireuit 
model  shown  in  Figure  8.  SF2.M  is  a  funetion  that  uses  the  differential  equations  (Equation  4) 
diseussed  previously  to  determine  the  state  variables.  Appendix  5  is  a  users  manual  and  explains 
in  more  detail  how  to  use  the  simulation. 
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Figure  24:  Bench  Generator  Simulation 


The  comparison  of  the  simulation  to  the  actual  lab  results  are  shown  in  Figure  25  as  well  as 
Table  3  below.  Also  shown  for  comparison  are  the  short  circuit  currents  as  calculated  with  the 
classical  equation.  The  results  match  within  0.1%  for  the  first  simulated  peak  compared  to  the 
actual  data.  The  second  short  circuit  peak  has  the  highest  error  and  then  decreases  to  -0.7%  at  42 
ms.  The  most  critical  bench  generator  parameter  to  match  is  the  first  short  circuit  peak 
amplitude  as  ultimately  the  interest  is  in  evaluating  short  circuit  protective  actions. 


51 


Figure  25:  Bench  Generator  (Simulation,  Classical,  and  Actual)  Results 


Iphase 

Peak  at  5  ms 

Peak  at  1 7  ms 

Peak  at  29  ms 

Peak  at  42  ms 

Actual  (A) 

170 

106 

72 

58 

Classical  (A) 

170.6 

0.4% 

112.0 

5.7% 

76.9 

6.8% 

57.0 

-1.7% 

Simulation  (A) 

169.8 

-0.1% 

109.9 

3.7% 

77.4 

7.5% 

57.6 

-0.7% 

Table  3:  Simulation  and  Classical  results  compared  to  Actual  Data 


5.1.2  Sensitivity  Analysis 

A  sensitivity  analysis  was  conducted  on  the  yellow  parameters  listed  in  Table  2  to  understand 
which  parameters  are  having  the  biggest  impact  on  the  equal  mutuals  model.  The  results  are 
shown  below  in  Table  4. 
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Sensitivity 

Actual 

Short  Circuit  Peak  Percent 

Short  Circuit  Peak  Percent 

Parameter 

Value 

Decrease  (%)  due  to  10%  decrease 

Increase  (%)  due  to  10%  increase 

in  sensitivity  parameter 

in  sensitivity  parameter 

Xq  (P.U.) 

0.77 

0.03 

0.05 

Xq  (p.U.) 

0.15 

1.3 

1.3 

x'd  (P-u.) 

0.197 

-0.5 

-0.6 

x'd  (p-u.) 

0.094 

-6.5 

-5.9 

tqo  (sec) 

0.07 

-0.05 

-0.03 

t'L  (sec) 

0.03 

0.7 

0.6 

Table  4:  Bench  Generator  Model  Sensitivity  Parameters  (First  Peak  Cnrrent) 


The  most  difficult  item  to  match  in  the  equal  mutuals  model  was  the  first  short  circuit  current 
peak  at  6  ms  corresponding  to  phase  b  current  as  shown  in  Figure  25.  This  was  a  peak  value  of 
170  amps  in  the  actual  short  circuit  test.  This  peak  value  for  phase  current  was  used  to  assess  the 
sensitivity  numbers  shown  in  the  table.  A  sensitivity  parameter  was  changed  either  plus  or 
minus  ten  percent  and  the  peak  current  percent  change  was  assessed.  A  negative  percentage 
change  corresponds  to  peak  current  increasing  for  a  decrease  in  a  sensitivity  parameter  or  vice 
versa.  The  most  sensitive  parameters  were  the  sub-transient  d-axis  and  q-axis  reactances 
respectively  (highlighted  in  yellow  in  the  table).  This  result  is  anticipated  based  on  the 
parameter  chosen  to  be  assessed  which  is  predominately  a  sub-transient  item  (i.e.  the  first  peak 
current). 

Another  sensitivity  analysis  was  conducted  on  the  decay  amplitude  at  42  ms.  This  was  looked  at 
to  determine  the  sensitivity  of  the  decay  envelope  of  the  short  circuit  current  as  a  result  of 
changing  the  same  parameters  listed  in  Table  2.  The  first  peak  amplitude  and  the  last  peak 
amplitude  for  the  short  circuit  currents  were  mutually  dependent  on  the  sensitivity  parameters 
identified. 
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Sensitivity 

Actual 

Short  Circuit  Envelope  Percent 

Short  Circuit  Envelope  Percent 

Parameter 

Value 

Decrease  due  to  10%  decrease  in 

Increase  due  to  10%  increase  in 

sensitivity  parameter 

sensitivity  parameter 

Xq  (P.U.) 

0.705 

-0.2 

-0.5 

Xq  (p.U.) 

0.243 

-0.4 

-0.1 

x'd  (P-u.) 

0.197 

-7.4 

-6.7 

x'd  (p-u.) 

0.0579 

3.7 

2.9 

tqo  (sec) 

0.0672 

0.9 

0.7 

t'L  (sec) 

0.048 

3.4 

3.5 

Table  5:  Bench  Generator  Sensitivity  Parameters  (Decay  Envelope) 


The  sub-transient  q-axis  open  circuit  time  constant  and  q-axis  reactance  had  little  effect  on  either 
the  first  peak  amplitude  for  short  circuit  current  or  the  decay  envelope  as  shown  in  Table  5.  The 
most  sensitive  parameter  for  affecting  both  the  first  peak  amplitude  and  the  decay  envelope  in 
the  simulation  was  the  sub-transient  d-axis  reactance  and  therefore  care  must  be  taken  when 
adjusting  this  parameter. 


5.1.2  DDG  Model 

The  equal  mutuals  model  developed  for  the  bench  generator  was  used  to  simulate  a  three  phase 
short  circuit  for  a  DDG  electrical  generator.  The  model  used  plant  parameters  for  two  DDG 
plants  as  provided  in  Table  6  below  to  run  two  different  simulations. 
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DDG-I 

DDG-II 

DDG-I 

DDG-II 

Phase  (KVAR) 

3125 

3125 

Xkdl(II)  p.U. 

0.333829 

0.400595 

^ base  (^)line-line  RMS 

450 

450 

Xq(a)  p.U. 

1.08 

1.296 

Ibase  (A)rmS 

4009 

4009 

Xd(Q)  p.U. 

1.848 

2.2176 

Pf 

0.8 

0.8 

Xkq(II)  p.U. 

1.3298 

1.5957 

Pole  Pairs 

2 

2 

Xq”(a)  p.U. 

0.328 

0.3936 

H  (sec) 

2.137 

1.7083 

Xd’(a)  p.U. 

0.207 

0.2484 

Ta  (Q)  p.U. 

0.00515 

0.00618 

Xd  (12)  p.U. 

0.172 

0.2064 

X,  (Q)  p.U. 

0.08 

0.096 

Xf(I2)  p.U. 

1.9048 

2.2858 

Xaq  (Q)  p.U. 

1 

1.2 

Xkd(I2)  p.U. 

2.1018 

2.5222 

Tkq  (Q)  p.U. 

0.0613 

0.07356 

tdo'  (sec) 

4.7846 

4.7846 

Xkql  (f^)  p.U. 

0.329787 

0.395744 

tqo”  (sec) 

0.10081 

0.10081 

Xad  (f2)  p.U. 

1.768 

2.1216 

tdo”  (sec) 

0.048521 

0.048521 

rf(a)  p.U. 

0.00111 

0.001332 

td  (sec) 

0.5553 

0.5553 

Xfl(a)  p.U. 

0.136829 

0.164195 

td  (sec) 

0.038912 

0.038911 

rkd(f2)  p.U. 

0.023968 

0.028762 

Table  6:  DDG-51  generator  parameters  provided  by  NAVSEA 


The  values  indicated  in  the  un-shaded  blocks  above  were  provided  by  NAVSEA  and  the  shaded 
parameters  were  calculated  from  the  provided  values  (Appendix  7).  The  equal  mutuals  model 
allows  for  the  generator  to  be  at  different  loading  levels  prior  to  the  short  circuit.  Two  separate 
simulations  were  run  with  the  DDG-I  parameters,  the  first  scenario  of  which  was  for  the 
generator  unloaded  and  the  second  scenario  was  for  the  generator  fully  loaded  at  rated  capacity. 
These  results  are  shown  in  Figure  26  and  Figure  27.  The  results  for  the  DDG-II  parameters  were 
similar  to  the  DDG-I  results  and  are  shown  in  Appendix  4. 
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Figure  26:  DDG-I  Short  Circuit  (luitially  Uuloaded) 


Figure  27:  DDG-I  Short  Circuit  (luitially  Fully  Loaded) 
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5.2  DDG  and  Bench  Generator  Comparison 


Referring  to  Figure  7  we  ean  observe  that  the  terminal  impedanee  will  have  a  seeond  order 
numerator  and  a  seeond  order  denominator.  The  terminal  impedance  can  be  written  in  terms  of 
the  transient  reactance  (xd),  transient  short  circuit  and  open  circuit  time  constants  (T^  and 
as  well  as  the  sub-transient  short  circuit  and  open  circuit  time  constants  (T^'  and  T^'g).  (6) 

(i+ns)a  +  7'M 

Equation  5:  Per  Unit  Inductance,  Synchronous  Generator 

From  section  4.1  it  was  observed  that  the  bench  generator  has  a  more  rapid  response  than  the 
DDG-51  generator.  The  bench  generator  transient  response  is  11.68-12.44  times  faster  than  the 
DDG  response  and  the  sub-transient  response  is  1.6  -  2.6  times  faster.  This  is  shown  in  the  pole 
zero  plot  in  Figure  28. 
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Figure  28:  Pole  Zero  Plots  for  Bench  Generator  and  DDG-I  Model 

From  these  results  and  the  impedance  magnitude  plots  in  Figure  29  it  is  possible  to  compare  the 
bench  generator  and  DDG  generator  responses  with  each  other. 
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Figure  29:  Impedance  Bode  Plots  (Bench  Generator  and  DDG-I) 


It  is  apparent  from  the  Bode  plots  that  the  DDG  generator  has  three  defined  regions  that  are 

Tad 

visible  in  Figure  29.  For  the  DDG  generator  up  to  approximately  0.2  is  the  steady  state 
region,  from  approximately  (3-10)  is  the  transient  region,  and  greater  than  approximately 

vcici 

60  is  the  sub-transient  regime.  The  beneh  generator  has  well  defined  steady  state  and  sub¬ 
transient  regions  but  there  is  not  a  well  defined  plateau  for  the  transient  regime.  This  is  not  a 
problem  for  the  modeling  being  accomplished  and  may  not,  in  fact,  be  a  problem  at  all.  If  the 
magnitude  response  is  determined  more  accurately  in  future  studies  by  conducting  a  frequency 
response  at  the  stator  terminals  it  would  be  possible  to  match  the  Bode  response  more  accurately 
in  the  future  by  adding  additional  poles  and  zeros  to  Equation  5.  For  the  purposes  of  the 
modeling  being  conducted  the  concern  is  not  to  match  the  transient  Bode  magnitude  response  but 
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instead  to  capture  the  sub-transient  response.  This  allows  us  to  look  at  fault  detection  times  and 
compare  those  to  fault  detection  times  for  the  DDG-5 1 . 


Figure  30:  DDG-51  Current  Threshold  after  Three  Phase  Fault  (Initially  Unloaded  -  Left),  (Initially  Fully 

Loaded  -  Right) 


Figure  31:  Bench  Generator  Current  Threshold  (Initially  Unloaded  -  Left),  (Initially  Fully  Loaded  -  Right) 
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Figure  30  and  Figure  31  show  the  length  of  time  for  a  three  phase  fault  to  cross  a  current 
threshold  of  150%  of  rated  capacity  for  both  the  bench  and  DDG  generators  (indicated  at  29.46 
and  8504  amps  respectively).  Both  of  these  generators  are  shown  with  initial  conditions  of  open 
circuited  terminals  (i.e.  unloaded)  as  well  as  the  case  of  a  generator  that  is  fully  loaded.  This 
150%  rated  current  threshold  was  used  as  a  reference  point  to  compare  the  DDG  and  bench 
generator  responses  to  the  approximate  amount  of  time  that  would  be  necessary  for  the  High 
Speed  Relay  (HSR)  algorithm  to  detect  a  fault  and  determine  whether  or  not  the  fault  was 
downstream  or  not.  This  was  used  in  lieu  of  a  power  calculation  since  the  phase  to  phase 
terminal  voltages  in  the  modeling  employed  was  zero  volts.  The  150%  current  threshold  is  a 
conservative  estimate  when  compared  to  the  power  threshold  used  in  the  HSR  algorithm  and  is  a 
good  datum  point  to  start  assessing  fault  detection  capabilities  for  NILM  testing.  The  figures 
above  show  all  three  phases  of  current  simultaneously  during  the  fault  scenarios.  The  times 
listed  for  the  fault  threshold  crossing  correspond  to  the  first  phase  of  current  which  crosses  this 
threshold.  The  bench  generator  and  DDG  currents  were  at  the  same  exact  phases  at  the  time  of 
the  simulated  three  phase  fault  and  therefore  can  be  related  to  each  other.  Two  observations  are 
that  number  one  the  threshold  is  crossed  faster  when  either  the  bench  or  DDG  generator  is  fully 
loaded  and  number  two  that  the  threshold  is  crossed  faster  for  the  bench  generator  compared  to 
the  DDG  generator  in  both  situations.  Figure  32  below  shows  actual  fault  data  corresponding  to 
the  bench  generator.  The  equal  mutuals  model  predicts  a  current  threshold  crossing  for  the 
bench  generator,  when  initially  unloaded,  of  0.69  milliseconds  whereas  the  actual  data  shows  a 
threshold  crossing  at  0.90  milliseconds.  This  difference  is  dependent  on  the  particular  phase  of 
current  that  is  captured.  Note  that  the  simulation  shows  each  phase  of  current  exceeding  the 
150%  limit  over  a  time  range. 
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Figure  32:  Bench  Generator  Actual  Phase  Current  Crossing  150%  Threshold 


Table  7  below  summarizes  the  times  to  exeeeding  the  eurrent  threshold  as  well  as  indieating  the 
initial  rate  of  eurrent  inerease.  For  the  beneh  generator  the  aetual  as  well  as  the  simulated  time  is 
shown  for  the  unloaded  eondition  sinee  that  data  point  was  available.  For  the  rest  of  the 
seenarios  just  the  simulated  time  is  shown. 


Time  (ms) 

Current  Rate  (A/ms) 

Bench  Generator  (Initially  Unloaded) 

0.90  (Actual) 

0.69  (Simulated) 

32.7  (Actual) 

42.7  (Simulated) 

Bench  Generator  (Initially  Fully  Loaded) 

0.21 

140.3 

DDG  (Initially  Unloaded) 

1.4 

6075 

DDG  (Initially  Fully  Loaded) 

0.41 

20742 

Table  7:  Generator  Times  for  Fault  Detection  (150%  Current) 
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5.2.1  Armature  Reaction  for  DDG  and  Bench  Generator 


The  armature  reaction  is  clearly  visible  for  the  case  of  the  DDG  simulation  but  decays  rapidly  for 
the  case  of  the  bench  generator  as  seen  in  Figure  24  thru  Figure  27.  The  equation  for  the 
armature  time  constant,  which  determines  how  rapidly  the  DC  offset  decays  to  zero,  is  shown 
below. 


T  = 

^  a 


h. 

(jOTn 


The  negative  sequence  reactance  (x2)  is  equal  to  *  X^.  The  armature  time  constants  for  the 
DDG  and  bench  generator  are  shown  in  Table  8  below. 


Ta  (H)  p.U. 

x'a(fi)  p.U. 

Xq  (Q)  p.U. 

X2  (D)  p.U. 

Ta  (ms) 

Bench  Generator 

.056 

.094 

.150 

.119 

5.6 

DDG 

.00515 

.172 

.328 

.238 

122.0 

Table  8:  Armature  time  constants  for  DDG  and  Bench  Generator 


To  illustrate  the  effects  of  the  armature  reaction  on  the  short  circuit  currents  the  simulation  in 
Appendix  4  was  run  on  the  DDG  parameters  for  three  cases.  The  first  case  used  the  actual  DDG 
parameters  while  the  second  and  third  case  used  armature  resistances  of  ten  times  and  twenty 
times  the  actual  DDG  armature  resistance  listed  in  Table  8.  This  had  the  effect  of  reducing  the 
armature  time  constant  to  12.2  and  6.1  milliseconds  respectively.  The  results  are  shown  in 
Figure  33  and  show  that  the  armature  reaction  is  another  important  difference  between  the  DDG 
and  bench  generators.  Whereas  the  bench  generator  short  circuits  are  primarily  a  result  of  the 
sub-transient  and  transient  reactances  and  time  constants  the  DDG  generator  also  has  a 
significant  armature  component. 
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Figure  33:  Effect  of  Armature  Resistauce  on  Short  Circuit  Current  for  DDG 


In  summary,  the  bench  generator  has  a  quicker  response  time  than  the  DDG  generator.  In  fact, 
for  the  loaded  condition  the  bench  generator  is  1.78  times  faster  than  the  DDG  whereas  for  the 
unloaded  case  the  bench  generator  is  1.53  times  faster.  This  means  that  when  fault  detection  and 
protective  action  testing  starts  to  be  evaluated  for  the  bench  generator  it  can  be  compared  to  the 
DDG  generator.  If  a  fault  detection  and  protective  action  can  be  shown  to  be  effective  for  the 
bench  generator  than  a  case  could  be  made  that  the  protective  action  would  be  effective  for  the 
DDG  plant  which  has  a  slower  response  time. 
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Chapter  6 


6.0  Future  Studies 

The  bench  generator  actual  lab  results  have  been  compared  to  the  equal  mutuals  model 
simulation  results  and  match  within  1.1%  for  the  three  phase  short  circuit  time  for  current  to 
cross  the  150%  threshold  set  point,  as  well  as  matching  initial  peak  current  amplitude  within 
0.1%.  Further  testing  and  validation  of  simulation  and  lab  results  are  required  to  validate  the 
bench  generator  model  but  the  results  so  far  are  encouraging.  Enough  validation  has  occurred  to 
this  point  start  to  start  building  and  testing  the  bench  generator  electrical  distribution  system. 
The  electrical  distribution  system  will  consist  of  three  bench  generators  that  will  supply  separate 
electrical  distribution  panels  that  are  separated  by  circuit  breakers.  The  setup  will  be  consistent 
with  a  typical  DDG  electric  plant  configuration. 

6.1  Additional  Work 

The  equal  mutuals  model  that  was  employed  consisted  of  one  set  of  damper  windings  for  the 
direct  and  quadrature  axes.  It  was  shown  in  chapter  four  that  the  direct  axis  terminal  did  not 
have  a  well  defined  transient  impedance  for  the  bench  generator  (Figure  29).  This  was  a  result  of 
the  transient  short  circuit  time  constant  being  relatively  close  to  the  sub-transient  open  circuit 
time  constant.  It  is  possible  to  apply  a  low  voltage  across  the  stator  terminals  and  vary  the 
frequency  from  10'  -  100  Hz  to  determine  the  frequency  response  of  the  synchronous  generator. 
This  can  be  done  with  the  rotor  aligned  to  both  the  direct  and  quadrature  axes  to  obtain 
information  about  both  Xq(s)  and  Xd(s)  (10).  Once  the  frequency  information  is  obtained  it 
would  then  be  possible  to  match  the  frequency  response  more  accurately  by  using  as  many 
damper  windings  as  required  as  opposed  to  the  single  winding  that  was  used  in  the  analysis 
presented  thus  far.  This  could  then  be  used  to  validate  the  results  presented  in  this  thesis.  It  is 
important  to  note  that  the  bench  generator  simulation  based  on  one  set  of  direct  and  quadrature 
windings  is  producing  simulation  results  that  match  the  sub-transient  response.  The  transient 
response  (i.e.  between  sub-transient  and  steady  state)  in  the  equal  mutuals  simulation  may  be 
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able  to  be  matehed  further  pending  the  results  from  the  frequeney  response  testing.  If,  however, 
the  results  of  the  frequency  response  testing  match  the  Figure  29Bode  results  for  the  bench 
generator  than  the  one  damper  winding  model  is  sufficient  for  the  simulation. 

The  next  important  step  is  to  use  the  electric  modeling  results  from  this  thesis  and  the  speed 
modeling  results  (12)  to  build  the  mock  DDG  electric  plant.  From  this  it  will  be  possible  to 
simulate  typical  loads  on  an  electrical  distribution  system  and  to  start  collecting  NILM  data  to  be 
used  for  fault  detection  analysis.  Finally,  a  system  of  circuit  breakers  set  up  to  simulate  a  DDG 
electric  plant  could  be  controlled  by  outputs  from  NILM  based  on  the  data  collected  from  the 
NILM  current  and  voltage  sensors.  The  NILM  outputs  would  be  programmed  to  trip  the  circuit 
breakers  to  protect  the  electrical  distribution  panels  based  on  manually  inserted  faults.  Data 
collected  would  be  used  to  determine  optimal  power  thresholds  to  determine  if  the  fault  is 
downstream  or  upstream. 

6.2  Conclusion 

The  synchronous  generator  modeling  accomplished  has  been  for  the  purpose  of  demonstrating 
that  testing  results  from  the  bench  generator  can  be  related  to  a  DDG  plant.  It  has  been  shown 
that  the  bench  generator  has  a  faster  electrical  response  than  the  DDG  plant  using  the  results 
from  the  equal  mutuals  model  simulations.  It  is  not  a  disadvantage  that  the  bench  generator  has 
faster  electrical  characteristics  than  the  DDG  plant  since  if  NILM  can  be  shown  to  provide  a 
benefit  for  fault  detection  and  protection  on  the  bench  generator  platform  it  stands  to  reason  that 
it  would  provide  the  same  benefit  to  the  slower  responding  DDG  plant.  A  main  benefit  of  the 
NILM  would  be  in  assisting  the  MFM  with  HSR  algorithm.  As  stated  previously  it  is  a  lemma 
of  this  thesis  that  NILM  could  be  used  to  provide  real  time  fault  threshold  tuning  for  a  DDG 
plant.  Chapter  five  results  have  shown  that  the  amount  of  time  to  cross  the  150%  rated  current 
capacity,  as  well  as  the  highest  current  peak  amplitude,  depends  on  the  initial  loading  of  the 
generators.  Currently,  MFM  with  HSR  algorithm  have  a  fault  detection  threshold  that  is  adjusted 
to  prevent  false  detections.  This  threshold  is  currently  set  independently  of  generator  loading, 
and  in  fact  must  be  set  for  the  worst  case  loading  condition  which  would  be  when  the  generators 
are  lightly  loaded.  Switching  transients,  which  occur  in  any  electrical  distribution  system,  are 
the  main  reason  why  fault  thresholds  need  to  be  set  higher.  NILM  could  be  used  to  identify 
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these  switching  transients  and  have  a  logic  system  that  prevents  them  from  being  classified  as 
faults.  This  would  allow  for  the  further  optimization  of  fault  set  points  by  not  having  to  set  the 
threshold  above  these  switching  transient  amplitudes.  The  next  stage  of  NILM  feasibility  for 
fault  detection  and  switchboard  protection  is  ready  to  commence  based  on  the  results  of  the 
generator  testing  and  simulations  discussed  in  this  thesis.  This  will  be  accomplished  by  testing 
on  a  mock  DDG  electric  plant  using  three  synchronous  bench  generators. 
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APPENDIX  1  -  MATLAB  Evaluation  of  Xd 


%  OCC  and  SCC  Curves  from  Raw  Lab  Data  called  with  OCCandSCC.m 
%  Gross  page  264  [7]  Electric  Machines 

%  This  matlab  code  does  not  call  any  other  functions  and  is  self  sufficient 
clear  all; 

%  Lab  Data 


%  OCCIf  and  OCCEa  are  laboratory  obtained  values 


OCCIf  =  [0  0.09  0.12  0.19  0.27  0.31  0.41  0.44  0.52  0.57  0.63  0.74  0.76  0.81  0.91  0.96  1  1.09 
1.14  1.24  ... 

1.24  1.32  1.39  1.46  1.52  1.59  1.63  1.68  1.75  1.79  1.87  1.92  1.96  2.03]; 

OCCEa  =  [0  21.4  25.04  34.7  45  53.1  67  72.3  84.1  91.6  100  112.9  115.9  122.5  133.7  139.2  142.7 
152.7  157.5  ... 

165.6  166.4  172.3  176.9  181.7  185.3  190.2  192  194.4  197.9  199.7  203.5  205.9  207.3  210.2]; 

%  SCCIf  and  SCCIa  are  laboratory  obtained  values 


SCCIf  =  [0  0.1  0.125  0.148  0.165  0.19  0.24  0.265  0.295  0.328  0.362  0.394  0.425  0.46  0.486  ... 
0.511  0.541  0.562  0.591  0.613  0.635  0.655  0.672  0.699  0.724  0.741]; 


SCCIa  =  [0  1.6  2  2.4  2.67  3.12  3.97  4.32  4.84  5.39  6.15  6.61  7.19  7.81  8.24  8.77  9.42  9.61  ... 
10.3  10.8  11.1  11.5  12  12.3  13.1  13.4]; 


%  Raw  Lab  Data  plot  OCC  and  SCC 
figure(l); 

[AX, HI, H2]=plotyy(OCCIf, OCCEa, SCCIf, SCCIa); 

%  The  following  code  just  labels  both  y  axes  and  sets  axis  limits 

set(get(AX(l),’Ylabel'),'Strmg',’E_{\phi\phi}  (V_{RMS})'); 
set(get(AX(2),'Ylabel'),'Strmg',’I_{\phm}  (A_{RMS})'); 


%  formatting  for  figure  1 

axes(AX(l)) 
axis([0  2.5  0  220]); 
grid; 

axes(AX(2)) 
axis([0  2.5  0  80]); 


%  handle  for  left  axis 
%  This  is  left  y  axis  range 

%  handle  for  right  axis 
%  This  is  left  y  axis  range 
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grid; 

xlabel(’I_f(A)'); 

%  title('OCC  and  SCC  characteristics'); 

%  OCC  fitted  eurve 

%  Obtain  a  linear  fit  curve  and  make  sure  starts  at  zero  therefore  mirror  about  (0,0) 

linearx  =  [-l*OCCIf(10:-l;l)  OCCIf(l:10)]; 
lineary  =  [-l*OCCEa(10;-l;l)  OCCEa(l:10)]; 

Kagtemp  =  polyfit(linearx,hneary,l); 

Kag  =  Kagtemp(l); 

xl  =  [0;.01;1.2];  %  Range  of  values 

y  1  =  xl *Kag;  %  Einear  y  values 

ind  =14;  %  ind  adjusted  until  fitted  curve  matehes 

last  =  length(OCCIf); 

Bx  =  log(  (OCCIf(last)-OCCEa(last)/Kag)/(OCCIf(ind)-OCCEa(ind)/Kag)  )/  (OCCEa(last)  - 
OCCEa(ind)); 

Ax  =  (  OCCIf(last)-  OCCEa(last)/Kag  )/exp(Bx*OCCEa(last)); 

Efl  =  0:. 1:210; 

Ifl  =  Efl/Kag  +  Ax*exp(Bx*Efl);  %  See  ehapter  4  of  thesis 


figure(2); 

plot(OCCIf,OCCEa,Ifl,Efl);  %  plot  laboratory  OCC  versus  fitted  OCC 

xlabel('I_f  (A)');  %  label  x  axis 

ylabel('E_{\phi\phi}  (V_{RMS})'); 

axis([0  2.03  0  210.2]); 

grid; 

%  title('I_f  =  E_{AB}/K_{ag}  +  A_xe^{B_xE_{AB}}'); 
text(1.4,80,'K_{ag}  =  165','EontSize',18); 
text(1.4,60,'A_x  =  .0023’,'EontSize’,18); 
text(1.4,40,'B_x  =  .0276','EontSize',18); 
legend('Eab  Data','Eitted  Equation'); 

%  SCC  fitted  eurve 

linearx  =  [-l*SCCIf  SCCIf(2:length(SCCIf))]; 
lineary  =  [-l*SCCIa  SCCIa(2:length(SCCIf))]; 

%  This  ensure  fitted  curve  starts  at  0,0 

KSCC  =  polyfit(hnearx,lineary,l); 


%  this  mirrors  short  eireuit  eurrent  about  0,0 
%  mirror  y  about  0,0 


%  linear  fit 
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Ifl  =  0;.001:.741; 

la  =  Ifl  *  KSCC(1,1);  %  Fitted  linear  curve  for  armature  current 

figure(3); 

plot(SCCIf,SCCIa,Ifl,Ia)  %  plot  actual  versus  fitted 
grid; 

%  Graph  formatting 

xlabel('I_f(A)'); 
ylabel(’I_{\phin}  (A_{RMS})'); 
legend('Lab  Data'/Linear  Fif); 

%  title('Short  Circuit  Characteristics  (Min  Dong)'); 
text(.l,12,'Slope  =  17.4  A/A','FontSize',18); 

%  Xd  curve  (Xd  corresponds  to  la  =  Irating  in  SCC 

Ef2  =  6:. 1:210.2; 

112  =  Ef2/Kag  +  Ax*exp(Bx*Ef2);  %  See  thesis  for  equation  explanation 
Ia2  =  If2*KSCC(l,l); 

Xd  =  Ef2./Ia2;  %  Here  is  the  direct  axis  reactance  over  entire  range  of  armature  current 

figure(4); 

plot(If2,Xd);  %  Graph  formatting 

grid; 

xlabel('I_f); 

ylabel('X_d'); 


Xdag  =  Kag/KSCC(l,l);  %  Air  gap  reactance,  this  is  in  the  linear  regions 

Irating  =  5000/(3*120);  %  Direct  axis  reactance  for  machine  is  defined  at  Current  rating  for 

%  machine 

If3  =  Irating/KSCC(l,l);  %  field  current  corresponding  to  armature  current  rating 

for  temp  =1:1  :length(If2)  %  Eind  index  of  field  current  corresponding  to  armature  rating 
if  If2(temp)<If3 
index  =  temp; 
end 
end 

Ef3  =  Ef2(index);  %  Armature  voltage  at  current  rating 

Xd  =  Ef3/Irating;  %  Direct  Axis  reactance  at  rated  armature  current 

Pb  =  5000;  %  base  power  (rating) 

Vb  =  120*sqrt(3);  %  base  voltage  (rating) 
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Zb  =  Vb^2/Pb; 


%  Base  resistance 


xd  =  Xd/Zb  %  per  unit  direct  axis  reactance 

%  title('X_d  Ming  Dong'); 
text(l,4,'X_d  =  8.70\Omega','FontSize',18); 
text(l,3,'X_{dag}  =  9.46\Omega','FontSize',18); 
text(l,2,'x_d  =  1.01  p.u.','FontSize',18); 
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APPENDIX  2  -  MATLAB  tachometer  to  frequency  code 


%  Speed  profile  for  short  eircuits 
% 

%  This  eode  requires  .MAT  files  for  inputs.  These  .MAT  are  obtained  from 
%  the  oseilloseope  data  eapture  eode.  The  .MAT  files  neeessary  are  from 
%  the  taehometer  output  wtiioli  in  our  eonfiguration  was  at  the  gear  to  the 
%  DC  motor.  The  output  is  the  speed  profile  of  the  generator. 

elear  all; 

load  JAN19_6.mat;  %  This  was  the  oseilloseope  data  eapture. 

%  Taehometer  data  are  voltages.  Does  not  matter  what  oseope 
%  settings  are.  Voltages  are  eonverted  in  this  eode  to  frequeney. 
takl  =  data(:,l);  %  The  first  eolumn  in  this  case  was  tachometer. 

%  the  variable  t  is  time  from  oseope. 

avgl  =  mean(takl);  %  This  average  will  be  used  to  compare  hall 

%  output  and  convert  to  0  or  1 . 

up=0;  %  up  and  down  are  counters. 

down=0; 

for  x=l:l:length(t) 
if  takl(x)  >  avgl 
takl(x)  =  1; 
up=up+l; 
else 

takl(x)=0; 
down=down+l; 
end 
end 

transl=0;  %  This  will  count  transitions  from  0  to  1  and  vice 

%  versa 

arryl  =  [];  %  This  will  contain  the  times  at  each  transition. 

%  As  well  as  Hi  or  Low. 

for  x=  1 : 1 ;  length(t)- 1 
if  takl(x)-takl(x+l)~=0 

arryl  =  [arryl ;t(x)  takl(x);t(x+l)  takl(x+l)]; 
transl  =  transl+1; 
end 
end 


%  Count  hi  voltages  and  low  voltages  and  convert  to  binary 
%  hi  voltage  converted  to  1 
%  low  voltages  converted  to  0 
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cycles  =  transl/32;  %  total  number  of  electrical  cycles  for  generator 
%  during  oscope  data  capture.  Just  informational. 


frql=[];  %  This  will  hold  mid  point  times  and  frequencies, 

for  x=2:4:length(arryl)-5 

frql  =  [frql;arryl(x,l)+((arryl(x+3,l)-arryl(x,l)))/2  l/12*(12/32)*... 

( 1 2/ 1 4)  *2/ (arry  1  (x+3 , 1  )-arry  1  (x,  1 ))] ; 
end 

pa  =  polyfit(frql(:,l),frql(:,2),12);  %  Fit  a  curve  to  frequency  data. 


ta=  .00001;.00001;.l; 

B  =  2*pi*polyval(pa,ta); 

%  This  plots  real  time  frequency  and  fitted  curve. 
figure(l) 

plot(frql(:,l),frql(:,2),ta,B/2/pi); 

grid; 

yhm([35  65]); 

save  demotorspeed.mat  frql;  %  Save  for  use  by  simulation. 
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APPENDIX  3  -  MATLAB  Bench  Generator  Model  Code 


%  Synchronous  Generator  Simulation 
%  Using  Notes  from  Professor  Kirtley's  Class  6.685 
%  Also  using  solution  set  for  PS  1 1  of  same  elass 
%  Uses  mi.m,  miq.m,  and  sf2.m  functions 

clear  all; 
elear  figures; 

global  ydd  ydk  ykd  ydf  yff  yqq  yqk  ykf  ykq  rf  rkd  rkq  omz  ra  rl  vf .. 
H  Tm  zx  endfreq  transientlength  etl  Ifprofile 


%%%%%%%%%%%%%%%%%%%%%%% 
%  INPUTS  % 

%%%%%%%%%%%%%%%%%%%%%%% 


%  Unknown  Values  -  Will  need  to  mateh  to  simulation 


xd=  1.0069; 
xq  =  .770486; 
xdp  =  .19676; 
xdpp  =  .094; 
xqpp  =.15; 
tdop  =  .4088; 
tdopp  =  .03; 
tqopp  =  .07; 
xl  = .05; 

H=  .19; 

%  Known  Values 


omz  =  60*2*pi; 
Pb  =  5000; 

Vb  =  120; 

Zb  =  Vb^2/Pb*3; 
Ib=Pb/Vb/3; 

%p=2; 

%rl=2592 


t0  =  0; 


%  synehronous  d-  axis  reaetanee 
%  synehronous  q-  axis  reaetanee 
%  transient  (d-axis)  reaetanee 
%  subtransient  d-  axis  reaetanee 
%  sub  transient  q-  axis  reaetanee 
%  transient  (open  eireuit)  time  constant 
%  subtransient  d-  axis  time  eonstant 
%  subtransient  q-  axis  time  eonstant 
%  armature  leakage  reaetanee 
%  rotor  inertia 


%  base  frequency 
%  base  power  (rating) 

%  base  voltage  (rating)  line  to  neutral 
%  We  are  208  line  to  line 
%  Impedanee  Base 
%  Current  Base 

%  Not  used  %  number  of  pole  pairs 
%  1296  Ohms  =  lOOW,  2592  Ohms  =  50W 
%  rl  used  if  not  doing  short  eireuit 

%  This  was  used  in  early  beneh  testing  with  resistor  values  used  in  lieu  of 
%  shorts 

%  Simulation  Start  Time  (seeonds) 
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tf  =  .08;  %  Simulation  Stop  Time 

[xad  xkdl  xfl  rkd  rf]  =  mi(xd,  xdp,  xdpp,  tdop,  tdopp,  xl,  omz);  %  Determine  direet  axis  model 

%  parameters 

Ifb=2;  %  Field  Current  Base 

Vfb=3*Vb*Ib/Ifb;  %  Field  Voltage  Base 
Zfb=Vfb/Ifb;  %  Field  Resistance  Base 

rf=15.8/Zfb;  %  Use  this  measured  resistance  ilo  calc  rf  in  mi.m  convert  to  per  unit 

xkd  =  xad  +  xkdl;  %  Based  on  circuit  model 
xf  =  xad  +  xfl;  %  Based  on  circuit  model 

xmd  =  [xd  xad  xad;  xad  xkd  xad;  xad  xad  xf]; 

ymd  =  inv(xmd); 

[xaq  xkql  rkq]  =  miq(xq,  xqpp,  tqopp,  xl,  omz); 
xkq  =  xaq  +  xkql; 
xmq  =  [xq  xaq;  xaq  xkq]; 
ymq  =  inv(xmq); 

ydd  =  ymd(  1,1);  %  inverted  xmd  matrix  parameters 

ydk  =  ymd(l,2); 
ykd  =  ymd(2,2); 
ydf=ymd(l,3); 
yff=ymd(3,3); 
ykf  =  ymd(2,3); 
yqq  =  ymq(l,l); 
yqk  =  ymq(l,2); 
ykq  =  ymq(2,2); 

ra  =  .48/Zb;  %  Measured  value  (use  this  one) 

vf  =  3 1/Vfb  %  Measured  vf  to  get  to  rated  voltage  (use  this  one) 

Tm  =  0;  %  Open  circuit  voltage  therefore  zero  Mechanical  Torque 

%  psiO  =  [psidO  psiqO  psikdO  psikqO  psifD  omO] 
psiO  =  [1010  1+xfl/xad  omz]; 

dt  =  1/8024;  %  Simulation  step  size 

time  =  t0;dt:tf;  %  Time  of  simulation 

load  ifO.mat;  %  This  is  real  field  current  from  short  circuit  test 

%  This  is  obtained  from  oscope. 

%  figure(l) 

%  plot(Ifprofde( : ,  1 )  *  1 000,lfprofile( :  ,2)); 


%  direct  axis  reactance  3x3  matrix  (equal 
%  mutuals  model) 

%  quadrature  axis  parameters 
%  q-axis  damping 

%  2x2  matrix  in  equal  mutuals  model 
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%  xlabel('Time  (ms)'); 

%  ylabel('Field  Current  (per  unit)'); 

%  grid; 

[t,psi]  =  ode45('sfIV[D',time,  psiO);  %  sfMD.m  fde 

id  =  ydd  .*  psi(:,l)  +  ydk  .*  psi(:,3)  +  ydf  .*  psi(:,5);  %  direct  axis  current 

iq  =  yqq  .*  psi(:,2)  +  yqk  .*  psi(:,4);  %  quadrature  axis  current 

Ib  =  Pb/(3*Vb);  %  Current  to  neutral  13.89  A 

a  =  2*pi/3; 

omz  =  psi(:,6);  %  omega  profile  factoring  in  inertia 

ia  =  sqrt(2)*Ib*(id  .*  cos(omz.*t)  -  iq  .*  sin(omz.*t)); 
ib  =  sqrt(2)*Ib*(id  .*  cos(omz.*t  -  a)  -  iq  .*  sin(omz.*t  -  a)); 
ic  =  sqrt(2)*Ib*(id  .*  cos(omz.*t  +  a)  -  iq  .*  sin(omz.*t  +  a)); 
iff  =  ydf.*  psi(:,l)  +  ykf  .*  psi(:,3)  +  yff  .*  psi(:,5); 
ikd  =  ydk  .*  psi(:,l)  +  ykd  .*  psi(:,2)  +  ykf.*  psi(:,3); 

figure(2); 

plot(t*  1 000,ia,t*  1 000,ib,t*  1 000, ic);  %  Line  to  neutral  currents 

title('Fault  Current,  SHORT  CIRCUIT  Simulation'); 

xlim([0  80]) 
ylim([-200  200]) 
ylabel('I_{\phi-n}  (A)'); 
xlabel('Time  (ms)'); 
legend('i_a','i_b','i_c'); 
grid 

MI.M 

%  model  elements  from  terminal  parameters  mi.m 

function  [xad,  xkd,  xf,  rkd,  rf]  =  mi(xd,  xdp,  xdpp,  tdop,  tdopp,  xl,  omz) 

xad  =  xd  -  xl; 

xf  =  xad  *  (xdp  -xl)  /  (xad  -  xdp  +  xl); 

xkd=  l/(l/(xdpp-xl)  -  l/xad-l/xf); 

rf  =  (xf  )/(omz  *  tdop); 

rkd  =  (xkd  +  xad*xf/(xad  +xf))/(omz*tdopp); 

MIO.M 

%  model  elements  from  terminal  parameters  miq.m 
function  [xaq,  xkq,  rkq]  =  miq(xq,  xqpp,  tqopp,  xl,  omz) 
xaq  =  xq  -xl; 
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xkq  =  xaq*(xqpp  -xl)/(xaq-xqpp+xl); 
rkq  =  (xaq+xkq)/(omz*tqopp); 


SF2.M 


function  dpsi  =  sf2(t,psi) 

global  ydd  ydk  ykd  ydf  yff  yqq  yqk  ykf  ykq  rf  rkd  rkq  omz  ra  rl  vf  H  Tm  zx... 
endfreq  transientlength  wO  Ifprofile  frql 

psid  =  psi(l); 
psiq  =  psi(2); 
psikd  =  psi(3); 
psikq  =  psi(4); 
psif  =  psi(5); 
om  =  psi(6); 

index  =  t/Ifprofde(end,l)*length(Ifprofde);  %  Used  to  index  to  aetual 

%  field  eurrent 

id  =  ydd*psid  +ydk*psikd  +  ydUpsif; 
ikd  =  ydk*  psid  +  ykd*psikd  +  ykf*psif; 

iff=Ifprofile(round(index)+l  ,2);  %  Aetual  Field  Current  Data 

iq  =  yqq  *  psiq  +  yqk*  psikq; 
ikq  =  yqk  *  psiq  +  ykq*  psikq; 


te  =  psid*iq  -  psiq  *id; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


psidotl  =  om*(psiq)  -  omz*(ra)*id  ; 
psidot2  =  -om*psid  -  omz*(ra)*iq; 
psidotS  =  -omz*rkd*ikd; 
psidotd  =  -omz*rkq*ikq; 
psidotS  =  omz*vf  -  omz*rf*iff; 

omdot  =  (omz/(2*H))*(te  +  Tm); 

dpsi  =  [psidotl  psidot2  psidotS  psidotd  psidotS  omdot]'; 
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PLOTSHORT.M 


%  plotshort.m 

%  This  m  file  determines  xdp  and  xdpp  by  plotting  elassical  equation  versus  %  %  oscilloseope 
output.  Also  determines  Tdp  and  Tdpp. 

elear  all 

load  w6inertia_JAN  1 9. mat;  %  radial  speed  profile  from  JAN  1 9  6. mat 
wreal=B;  %  B  is  from  w6inertia_JAN19.mat 

%  This  is  used  in  elassieal  equation 

load  JAN19_6.mat;  %  ehannel  1  is  taehometer 

%  Channel  2  is  Vswiteh 
%  Channel  3  is  Field  Current 
%  Channel  4  is  Short  eireuit  Current 


V  =  120*sqrt(2);  %  used  for  elassieal  plot  Voltage  to  neutral 

delay  =  pi/1 80*200;  %  Mateh  elassical  phase  to  oscope  and 

%  simulation.  Purely  to  make  sine  waves  be  on 
%  same  phase  (visually). 

Xd  =  8.7;  %  Known  value.  See  Appendix  I. 

%  Xdp,  Xdpp,  Tdp,  and  Tdpp  are  estimated  from  this  program  (trial  and 
%  error) 

Xdp  =1.7;  %  Match  Xdp,  Xdpp,  Tdp,  and  Tdpp  to  oscope  data 

Xdpp  =  .81216; 

Tdp  =  45e-3; 

Tdpp  =15e-3; 

ib  =  [];  %  Capture  classical  equation  current 

t0=18.3e-3;  %  When  current  starts  (match  delay  from  oscope  to 

%  classical  equation).  When  switch  turned  on. 

%  to  is  based  on  looking  at  oscilloscope  real  data.  When  does  short  start. 

for  x=l:l:length(t) 
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ift(x)  <  to  %  Ib  ILO  Vswitch  activated  since  not  at  same  time 

ib  =  [ib;0]; 
else 


ib  =  [ib;V*(l/Xd  +  (1/Xdp  -l/Xd)*exp(-(t(x)-tO)/Tdp)  +  (1/Xdpp  -  l/Xdp)*exp(-(t(x)- 
toyidpp))... 

*sin(wreal(x)  *  t(x)-delay)] ; 
end 
end 


figure(3); 

plot((t)*1000-18.3,data(:,4)*1000,(t)*1000-18.3,ib,simdata(:,l),simdata(:,2),t*1000,- 
ones(length(t),  l)*threshold); 
ylim([-200  200]); 

grid; 

xlim([0  50]); 
xlabel('Time  (ms)'); 
ylabel('I_{\phin}  (A)'); 

legend('I_{\phin}  FLUKE', 'I_{\phin}  Classical', 'I_{\phin}  Simulation'); 
title('Short  Circuit  Current'); 
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APPENDIX  4  -  MATLAB  DDG-51  Simulation  Code  and  DDG-II  Results 


%  DDG  Simulation 

%  This  code  uses  inputs  from  NAVSEA  to  simulate  a  three  phase  symmetrical  %  %  %  short 
clear  all; 

global  ydd  r  1  ydk  ykd  ydf  yff  yqq  yqk  ykf  ykq  rf  rkd  rkq  omz  ra  vf  Tm  H  delta 
%  DDG  generator  1 
Pb  =  3125e3; 

Vb  =  450/sqrt(3);  %  line  to  neutral 

Zb  =  Vb^2/Pb*3;  % 

Ib=PbA^b/3; 
pf=  .8; 

p  =  2; 

H  =  2.137; 
ra  =  .00515; 
xl=  .08; 

xaq  =  1 ; 
rkq  =  .0613; 
xkql  =  .329787; 

xad  =  1.768; 
rf=  .00111; 
xfl  =  .  136829; 
rkd  =  .023968; 
xkdl  =  .333829; 

%  Calculations  to  determine  other  parameters 

xq  =  xaq  +  xl; 
xd  =  xad  +  xl; 

xkq  =  xaq  +  xkql; 

xqpp  =  xl  +  (xaq*xkql)/(xkql  +  xaq); 
xdp  =  xl  +  (xad*xfl)/(xad  +  xfl); 
xdpp  =  l/(l/xad  +  1/xkdl  +  l/xfl)+  xl; 
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omz  =  2*pi*60; 
xf  =  xad  +  xfl; 
xkd  =  xad  +  xkdl; 

tdop  =  l/(omz*rf)*(xfl  +  xad)+l /(omz  *rkd)*  (xkdl  +  xad); 

tqopp  =  l/(omz*rkq)*(xkq  +  xaq); 

tdopp  =  1 /(omz  *rkd)*  (xkdl  +  (xad*xfl)/(xad  +  xfl))... 

/(I  +(l/(omz*rkd)*(xkdl  +  (xad)))/(l/(omz*rf)*(xfl  +  (xad)))); 

tdp  =  l/(omz*rf)*(xfl  +  (xad*xl)/(xad  +  xl))  +  l/(omz*rkd)*(xkdl  +  (xad*xl)/(xad  +  xl)); 

tdpp  =  1 /(omz  *rkd)*  (xkdl  +  (xad*xfl*xl)/(xad*xfl  +  xad*xl  +  xfl*xl))... 

/(I  +(l/(omz*rkd)*(xkdl  +  (xad*xl)/(xad  +  xl)))/(l/(omz*rf)*(xfl  +  (xad*xl)/(xad  +  xl)))); 


%  Initial  conditions 

vt=l; 

%  terminal  voltage 

ia  =  0; 

%  terminal  eurrent  magnitude 

psil  =  acos(.8); 

%  power  faetor  angle 

i_a  =  ia*(pf  -  j*sin(psil)); 

%  as  a  eomplex  veetor 

V  i  =  vt  +  ra*i  a; 

%  this  is  the  internal  voltage  veetor 

vi  =  abs(v_i); 

%  and  this  is  the  absolute  value  of 
%  internal  voltage 

ai  =  angle(v_i); 

%  and  this  is  the  angle  of  the  internal 
%  voltage 

e_l  =  v_i  +  j*i_a*xq; 
el  =  abs(e_l); 

%  the  voltage  that  defines  the  q  axis 

delta  =  angle(e_l); 

%  and  this  is  the  phase  angle 

id  =  ia*sin(  delta  +  psil); 

%  this  is  direet  axis  eurrent 

iq  =  ia*eos(delta  +  psil); 

%  and  this  is  quadrature  axis  eurrent 

vq  =  vi*eos(delta  -  ai); 

%  q  axis  eomponent  of  that  internal 
%  voltage 

vd  =  vi*sin( delta  -  ai); 

%  d  axis  eomponent 

eaf  =  vq  +  id*xd; 

%  voltage  behind  synehronous  reaetanee 

i_f  =  eaf/xad; 

%  steady  field  eurrent 

vf  =  i  f*rf; 

%  exeiter  voltage 

psidO  =  vq; 

%  fluxes  are  -90  degrees  rotated  from 
%  internal  voltage 

psiqO  =  -vd; 

psikdO  =  psidO  +  id*xl; 

%  same  as  flux  on  the  magnetizing 
%  induetanee 

psikqO  =  psiqO  +  iq*xl; 

%  this  is  the  same  too 

psifD  =  psikdO  +  xfl*i_f; 

%  flux  initial  eondition 

Tm  =  psidO *iq  -  psiqO* id; 

%  generating  torque 
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psia  =  [psidO  psiqO  psikdO  psikqO  psifD  omz]; 
xmd  =  [xd  xad  xad;  xad  xkd  xad;  xad  xad  xf]; 
ymd  =  inv(xmd); 

xmq  =  [xq  xaq;  xaq  xkq]; 

ymq  =  inv(xmq); 

ydd  =  ymd(l,l); 

ydk  =  ymd(l,2); 

ykd  =  ymd(2,2); 

ydf=ymd(l,3); 

yff=ymd(3,3); 

ykf  =  ymd(2,3); 

yqq  =  ymq(l,l); 

yqk  =  ymq(l,2); 

ykq  =  ymq(2,2); 

%  END  INITIAL  CONDITIONS 

time  =  -.05: 1/2024;- 1/2024;  %  negative  values  used  to  annotate  pre-fault 

[tp,psip]  =  ode45('sfr,time,  psia); 

idp  =  ydd  .*  psip(:,l)  +  ydk  .*  psip(:,3)  +  ydf  .*  psip(:,5); 
iqp  =  yqq  .*  psip(:,2)  +  yqk  .*  psip(:,4); 


%  psiO  =  [psidO  psiqO  psikdO  psikqO  psifO  omzO] 

psiO  =  [psip(end,l)  psip(end,2)  psip(end,3)  psip(end,4)  psip(end,5)  omz]; 

t0  =  0; 
tf=.50; 


%  psiO  =  [psidO  psiqO  psikdO  psikqO  psifO  omz]; 
dt  =  (tf-t0)/2024; 

time  =  tO;dt:tf; 

[t,psi]  =  ode45('sf2',time,  psiO); 

id  =  ydd  .*  psi(:,l)  +  ydk  .*  psi(:,3)  +  ydf  .*  psi(:,5); 
iq  =  yqq  .*  psi(:,2)  +  yqk  .*  psi(:,4); 
tel  =  (psi(:,l).*iq  -  psi(:,2).*id); 

id  =  [idp;  id]; 
iq  =  [iqp;iq]; 
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t=  [tp;t]; 


a  =  2*pi/3; 

oma  =  [psip(:,6);psi(:,6)]; 

ia  =  id  .*  cos(oma.*t)  -  iq  .*  sin(oma.*t); 

ib  =  id  .*  cos(oma.*t  -  a)  -  iq  .*  sin(oma.*t  -  a); 

ic  =  id  .*  cos(oma.*t  +  a)  -  iq  .*  sin(oma.*t  +  a); 

iff  =  ydf  .*  psi(:,l)  +  ykf  .*  psi(:,3)  +  yff  .*  psi(:,5); 

ikd  =  ydk  .*  psi(:,l)  +  ykd  .*  psi(:,2)  +  ykf.*  psi(:,3); 

threshold  =  Ib*sqrt(2)*ones(length(t),l); 

figure(l) 

clf 

Ib  =  Pb/(3*Vb); 

plot(t*1000,  ia*Ib/1000*sqrt(2), 

ib*Ib/1000*sqrt(2), t*  1000, ic*Ib/1000*sqrt(2),t*  1000, threshold*  1.5/1000 
threshold*  1.5/1 000);%,  t,  ib,  t,  ic); 

%  title('Fault  Current  DDG'); 
ylabel('I_{line}  (KA)'); 
xlabel('Time,  ms'); 
xhm([-50  500]); 
grid  on; 

legend('I_A','I_B','I_C'); 

imax=max(ic  *  Ib) ; 

%  figure(2); 

%  subplot  211 
%  plot(t,  iq*Ib);grid; 

%  title('I_q  and  I_d'); 

%  subplot  212 
%  plot(t,id*Ib);grid; 

%  xlabel('Time  (seeonds)'); 


figure(3) 

plot(t*1000,oma/omz*60);grid; 
xlabel('Time  (ms)'); 
ylabel('Frequenoy  (Hz)'); 

%  Ifb  =  1.5*Vb*Ib/vf; 

%  figure(4) 

%  plot(t,  iff*Ib);grid; 

% 


t*1000, 

,t*1000,- 
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%  figure(5); 

%  Tb  =  .0658e6*omz/(2*H); 
%  plot(t,tel*Tb);grid; 

% 

%  figure(6); 

%  plot(t,psi(:,7)/pi*180);grid; 


DDG  II  Simulation  -  These  are  simulation  results  using  DDG-II  parameters  from  Chapter  5 
(Note  that  eode  above  shows  DDG-I  parameters) 


Figure  34:  DDG-II  Short  Circuit  Simulatiou  (luitially  Uuloaded) 
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60 


Time  (sec) 

Figure  35:  DDG-II  Short  Circuit  Simulatiou  Frequeucy  Respouse  (Initially  Unloaded) 
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APPENDIX  5  -  SIMULATION  USERS  GUIDE 


All  data  files  and  simulation  code  is  located  on  the  bucket  account  under  leghorn. 

Bench  Generator  Simulation  Notes 


PLOTSHORT.M 

1.  Plotshort.m  file  is  used  to  determine  the  short  circuit  transient  and  sub-transient  reactances 
and  time  constants.  X^,  T^,  T^) 

2.  These  inputs  are  entered  based  on  some  assumptions.  The  assumptions  for  the  generator 
being  tested  were  that  the  transient  reaction  was  dominant  compared  to  sub-transient 
reaction.  Therefore  the  initial  matching  technique  was  to  assume  zero  sub-transient  reaction 
and  match  transient  parameters  as  best  as  possible.  Next  match  sub-transient  parameters. 

3.  The  speed  profile  used  in  the  classical  equation  was  obtained  from  the  tachometer  results 
which  were  a  .MAT  file  from  Appendix  2  code.  Alternatively  the  speed  profile  from  the 
simulation  can  be  used  since  this  is  ultimately  matched  to  the  tachometer  frequency. 

MDSIMULATION.M 

1 .  The  functions  used  by  this  simulation  are  MI.M,  MIQ.M,  and  SFMD.M. 

2.  MI.M  calculates  direct  axis  parameters. 

3.  MIQ.M  calculates  quadrature  axis  parameters. 

4.  SFMD.M  is  the  function  used  by  ode45  and  is  based  on  the  equal  mutuals  model. 

5.  SFMD.M  uses  actual  data  for  field  current.  This  is  what  the  Ifprofile  parameter  is.  Field 
current  is  calculated  for  the  DDG  simulations. 

6.  The  rotor  inertia  constant  (H)  is  determined  by  matching  the  simulation  frequency  response 
to  the  actual  generator  slow  down.  The  larger  H  is  the  less  slowdown  will  be  seen  in  the 
simulation. 
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7.  The  armature  leakage  reactance  is  somewhat  arbitrary  but  it  is  important  to  make  sure  that 
none  of  the  equal  mutuals  parameters  are  negative  as  a  result  of  this  value.  Typical  range  of 
leakage  reactance  is  0.04  -  0.15. 

8.  There  are  typical  ranges  of  values  for  and  tqg)  (11)(9).  These  are  used  in 

conjunction  with  the  values  obtained  in  plotshort.m  for  an  initial  estimate.  The  sensitivity 
analysis  conducted  in  chapter  5  should  be  looked  at  to  determine  the  parameters  that  need  to 
be  adjusted  carefully. 

9.  The  simulation  can  now  be  adjusted  to  match  the  actual  oscilloscope  results. 

10.  IfO.mat  is  the  actual  field  current  data  from  the  oscilloscope. 

1 1 .  Initial  conditions  can  be  entered  by  adjusting  ia.  For  open  circuited  initially  this  is  set  to  zero 
and  for  fully  loaded  this  is  set  to  one.  Any  other  loading  condition  can  be  entered  also. 

12.  The  DDG  simulation  does  not  require  MI.M  or  MIQ.M  since  all  the  parameters  can  be 
directly  calculated  without  those  functions. 
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APPENDIX  6  -  MATLAB  MODEL  VERIFICATION 


The  MATLAB  model  employed  for  this  thesis  was  used  to  run  the  simulation  of  a  three  phase 
fault  using  the  parameters  shown  in  Figure  36.  The  actual  results  are  shown  in  Figure  37.  The 
simulation  model  employed  and  the  results  from  that  simulation  are  shown  below.  There  are  no 
specific  values  to  be  compared  but  the  envelopes  from  each  of  the  graphs  (Figure  38  -  Figure  40) 
match  the  Krause  results.  Also  the  magnitudes  match.  Based  on  the  matching  results  the  equal 
mutuals  simulation  model  is  validated.  It  should  be  noted  that  since  there  were  two  quadrature 
axis  damper  windings  provided  in  the  data  in  Figure  36  the  equal  mutuals  model  was  modified  to 
account  for  the  second  damper  winding.  This  also  shows  that  if  the  bench  generator  is  modeled 
with  more  than  one  damper  winding  this  can  easily  be  adapted  into  the  simulation  code. 


Table  5,10-2  Slcam  turbine  generator 

Riiiin{>;  S35  MVA 
(.inc  (inc  voltage;  26  kV 
Power  faelor;  ().S5 
Polc.s;  2 

Spccil:  .1600  r  min 

Combined  inertia  of  generator  and  turbine: 

./  0.065s  X  lO”  J  s",  or  It  /<’  -  1.56  x  10'’  Ibm  ft’  II  5.6  s 

f’arametcrs  in  ohms  and  per  unit: 

0.00.'’4.1  £2.  't.tXt.l  pu 
.V,,  0.1518  £2.  0.19  pu 

.V,,  -  1.457  £2.  1.8  pu  .Y,,  1.457  £2.  1.8  pu 

''In  ^  <>  £’‘044  £2.  O.tKil  78  pu  r',,,  =  0.(K)075  £2.  0.000929  pu 

•Luo  •  0.6578  £2.0.8125  pu  .V;,^  O.  M4  5  £2.  0.  (4 14  pu 

().t)0681  £2.  0.00841  pu  0.01080  £2.  0.01 3.14  pu 

•Yr,:  -  0.07602  £2.  0.0939  pu  .\  r.,  0.06577  £2.  0.08 125  pu 


Figure  36:  Krause  Simulatiou  used  to  verify  Equal  Mutuals  Simulatiou  Model 
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Figure  37:  Krause  Simulatiou  Results 


%  Krause  Simulation  from  Chapter  5 

%  This  was  used  to  validate  our  model.  Known  inputs  were  entered  into  the 
%  simulation  and  the  results  matched  known  outputs.  These  were  provided 
%  from  "Analysis  of  Electric  Machinery".  See  Bibliography. 

clear  all; 

global  ydd  ydk  ykd  ydf  yff  yqq  yqkl  yqk2  ykf  ykql  ykq2  y3  rf  rkd  rkql  rkq2  omz  ra  vf  Tm  H 
ChapSparams; 

ta  =  .5*(xdpp+xqpp)/(omz*ra); 

t0  =  0; 
tf=.362; 

xmd  =  [xd  xad  xad;  xad  xkd  xad;  xad  xad  xf]; 
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ymd  =  inv(xmd); 


xmq  =  [xq  xaq  xaq;  xaq  xkql  xaq;xaq  xaq  xkq2]; 

ymq  =  inv(xmq); 

ydd  =  ymd(l,l); 

ydk  =  ymd(l,2); 

ykd  =  ymd(2,2); 

ydf=ymd(l,3); 

yff=ymd(3,3); 

ykf  =  ymd(2,3); 

yqq  =  ymq(l,l); 
yqkl  =ymq(l,2); 
yqk2  =  ymq(l,3); 
ykql  =  ymq(2,2); 
ykq2  =  ymq(3,3); 
y3  =  ymq(2,3); 

%  flux  vector  is  [psid  psiq  psikd  psikq  psif] 

psiO  =  [psidO  psiqO  psikdO  psikq  1  psifD  omz  psil  psikq2]; 

dt  =  (tf-t0)/1024; 

Tm  =.852;  %  added  from  chapt  5 

time  =  tO;dt:tf; 

[t,psi]  =  ode45('sf  ,time,  psiO); 

id  =  ydd  .*  psi(:,l)  +  ydk  .*  psi(:,3)  +  ydf  .*  psi(:,5); 
iq  =  yqq  .*  psi(:,2)  +  yqkl  .*  psi(:,4)  +  yqk2*psi(:,8); 

tel  =  (psi(:,l).*iq  -  psi(:,2).*id); 


a  =  2*pi/3; 
oma  =  psi(:,6); 

ia  =  id  .*  cos(oma.*t)  -  iq  .*  sin(oma.*t); 
ib  =  id  .*  cos(oma.*t  -  a)  -  iq  .*  sin(oma.*t  -  a); 
ic  =  id  .*  cos(oma.*t  +  a)  -  iq  .*  sin(oma.*t  +  a); 
iff  =  ydf.*  psi(:,l)  +  ykf.*  psi(:,3)  +  yff  .*  psi(:,5); 
ikd  =  ydk  .*  psi(:,l)  +  ykd  .*  psi(:,2)  +  ykf.*  psi(:,3); 

%  open  circuit  out  of  phase 


psiO  =  psi; 

%  psiO  =  [00  00  1+xfl/xad]; 
tOl  =  .36201; 


91 


tfl  =  2.4; 

dt  =  (tfl-t01)/1024; 

timel  =  t01:dt:tfl; 
psi(end,7) 

[ta,psia]  =  ode45('sfr, timel,  psi(end,:)); 

ida  =  ydd  .*  psia(:,l)  +  ydk  .*  psia(:,3)  +  ydf  .*  psia(:,5); 
iqa  =  yqq  .*  psia(:,2)  +  yqkl  .*  psia(:,4)+  yqk2*psi(:,8); 

te2  =  -(psia(:,l).*iqa  -  psia(:,2).*ida); 


a  =  2*pi/3; 
oma  =  psi(:,6); 

ial  =  ida  .*  cos(oma.*ta)  -  iqa  .*  sin(oma.*ta); 
ibl  =  ida  .*  cos(oma.*ta  -  a)  -  iqa  .*  sin(oma.*ta  -  a); 
icl  =  ida  .*  cos(oma.*ta  +  a)  -  iqa  .*  sin(oma.*ta  +  a); 
iffl  =  ydf.*  psia(:,l)  +  ykf  .*  psia(:,3)  +  yff  .*  psia(:,5); 
ikdl  =  ydk  .*  psia(:,l)  +  ykd  .*  psia(:,2)  +  ykf.*  psia(:,3); 

%  group  terms 


ia  =  [ia;  ial]; 
ib  =  [ib;  ibl]; 
ie  =  [ie;  iel]; 
iff=  [iff;  iffl]; 
ikd  =  [ikd;  ikdl]; 
t  =  [t;  ta]; 
tel  =  [tel;  te2]; 
psi  =  [psi;  psia]; 
id  =  [id;  ida]; 
iq  =  [iq;  iqa]; 


figure(l) 

elf 

subplot  311 
lb  =  Pb/(sqrt(3)*Vb); 
plot(t,  -ib*lb);%,  t,  ib,  t,  ie); 
title('Fault  Current'); 
ylabel('l_{\phi}  (A)'); 

grid  on; 

imax=max(io*lb) 
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subplot  312 
plot(t,  -iq*Ib);grid; 
ylabel('I_q(A)'); 
subplot  313 
plot(t,-id*Ib);grid; 
xlabel('Time,  s'); 
ylabel('I_d  (A)'); 


figure(2) 
subplot  211 

plot(t,psi(:,6)/omz);grid; 
ylabel('\omega  (per  unit)'); 
subplot  212 

plot(t,psi(:,7)/pi*180);grid; 
ylabel('\theta  (degrees)'); 
xlabel('Time  (sec)'); 

figure(3) 
subplot  211 
plot(t,  iff*Ib);grid; 
ylabel('Field  Current  (A)'); 

subplot  212 

Tb  =  .0658e6*omz/(2*H); 
plot(t,tel  *Tb);grid; 
ylabel('Torque  (Nm)'); 
xlabel('Time  (sec)'); 


CHAP5PARAMS.M 


%  Parameters  for  Chapter  5  Simulation; 

Pb  =  835e6; 

Vb  =  26e3; 

Zb  =  Vb^2/Pb; 
pf  =  0.85; 

p  =  1 ;  %  Number  of  pole  pairs 

H  =  5.6; 

ra  =  .003; 
xl=  .19; 
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xq  =  1.8; 
rkql  =  .00178; 
xlkql  =  .8125; 
rkq2  =  . 00841; 
xlkq2  =  .0939; 


xkql  =  imag((rkql  +  j*(xlkql))*(rkq2  +  j*(xlkq2))/((rkql  +  j*(xlkql))  +  (rkq2  +  j*(xlkq2)))); 
rkq  =  real((rkql  +  j*(xlkql))*(rkq2  +  j*(xlkq2))/((rkql  +  j*(xlkql))  +  (rkq2  +  j*(xlkq2)))); 

xd  =  1.8; 
rf  =  .000929; 
xfl  =  .1414; 
rkd  =  .01334; 
xkdl=  .08125; 

xaq  =  xq  -  xl; 
xad  =xd  -  xl; 

xkq  =  xaq  +  xkql; 

xkql  =  xaq  +  xlkql;  %  Added 

xkq2  =  xaq  +  xlkq2;  %  Added 


xqp  =  xl  +  (xaq*xlkql)/(xlkql  +  xaq); 
xdp  =  xl  +  (xad*xfl)/(xad  +  xfl); 

xqpp  =  xl  +  (xaq*xlkql*xlkq2)/(xaq*xlkql  +  xaq*xlkq2  +  xlkql  *xlkq2); 
xdpp  =  xl  +  (xad*xfl*xkdl)/(xad*xfl  +  xad*xkdl  +  xfl*xkdl); 

omz  =  2*pi*60; 


tdop  =  l/(omz*rf)*(xfl  +  xad)+l/(omz*rkd)*(xkdl  +  xad); 

tqopp  =  l/(omz*rkq2)*(xlkq2  +  (xaq*xlkql)/(xlkql  +  xaq))... 

/(I  +(l/(omz*rkq2)*(xlkq2  +  xaq))/(l/(omz*rkql)*(xlkql  +  (xaq)))); 
tdopp  =  l/(omz*rkd)*(xkdl  +  (xad*xfl)/(xad  +  xfl))... 

/(I  +(l/(omz*rkd)*(xkdl  +  (xad)))/(l/(omz*rf)*(xfl  +  (xad)))); 


tqp  =  l/(omz*rkql)*(xlkql  +  (xaq*xl)/(xl  +  xaq))+l/(omz*rkq2)*(xlkq2  +  (xaq*xl)/(xl  +  xaq)); 
tdp  =  l/(omz*rf)*(xfl  +  (xad*xl)/(xad  +  xl))  +  l/(omz*rkd)*(xkdl  +  (xad*xl)/(xad  +  xl)); 
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tqpp  =  l/(omz*rkq2)*(xlkq2  +  (xaq*xlkql*xl)/(xaq*xl  +  xaq*xlkql  +  xl*xlkql))... 

/(I  +(l/(omz*rkq2)*(xlkq2  +  (xaq*xl)/(xaq  +  xl)))/(l/(omz*rkql)*(xlkql  +  (xaq*xl)/(xaq  + 

xl)))); 


tdpp  =  l/(omz*rkd)*(xkdl  +  (xad*xfl*xl)/(xad*xfl  +  xad*xl  +  xfl*xl))... 

/(I  +(l/(omz*rkd)*(xkdl  +  (xad*xl)/(xad  +  xl)))/(l/(omz*rf)*(xfl  +  (xad*xl)/(xad  +  xl)))); 
xf=xad  +  xfl;  % 

xkd  =  xad  +  xkdl; 


%  Initial  conditions 

vt=l; 
la  =  1; 

psil  =  acos(pf); 

i_a  =  ia*(pf  -  j*sin(psil)); 

v_i  =  vt  +  ra*i_a; 

vi  =  abs(v_i); 

ai  =  angle(v_i); 

e_l  =  v_i  +  j*i_a*xq; 

el  =  abs(e_l); 

delta  =  angle(e_l); 

id  =  ia*sin(  delta  +  psil); 

iq  =  ia*eos(delta  +  psil); 

vq  =  vi*eos(delta  -  ai); 

vd  =  vi*sin( delta  -  ai); 

eaf  =  vq  +  id*xd; 

i_f  =  eaf/xad; 

vf  = 

psidO  =  vq; 

psiqO  =  -vd; 

psikdO  =  psidO  +  id*xl; 

psikql  =  psiqO  +  iq*xl; 

psikq2  =  psiqO  +  iq*xl; 

psifD  =  psikdO  +  xfl*i_f; 


%  terminal  voltage 
%  terminal  current  magnitude 
%  power  factor  angle 
%  as  a  eomplex  veetor 
%  this  is  the  internal  voltage  veetor 
%  and  this  is  the  absolute  value  of  internal  voltage 
%  and  this  is  the  angle  of  the  internal  voltage 
%  the  voltage  that  defines  the  q  axis 

%  and  this  is  the  phase  angle 
%  this  is  direet  axis  eurrent 
%  and  this  is  quadrature  axis  eurrent 
%  q  axis  eomponent  of  that  internal  voltage 
%  d  axis  eomponent 

%  voltage  behind  synehronous  reaetanee 
%  steady  field  eurrent 
%  exeiter  voltage 

%  fluxes  are  -90  degrees  rotated  from  internal  voltage 

%  same  as  flux  on  the  magnetizing  induetanee 
%  this  is  the  same  too 

%  flux  initial  eondition 


Tm  =  psid0*iq  -  psiq0*id  %  generating  torque 


SF.M 

funetion  dpsi  =  sf(t,psi) 

global  ydd  ydk  ykd  ydf  yff  yqq  yqkl  yqk2  ykf  ykql  ykq2  y3  rf  rkd  rkql  rkq2  omz  ra  vf  H  Tm% 
xd  xf  xkd  xad  xaq  xq  xkq  ymd 


psid  =  psi(l); 
psiq  =  psi(2); 
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psikd  =  psi(3); 
psikql  =  psi(4); 
psif  =  psi(5); 
om  =  psi(6); 
delta  =  psi(7); 
psikq2  =  psi(8); 


id  =  ydd*psid  +ydk*psikd  +  ydf'psif; 
ikd  =  ydk*  psid  +  ykd*psikd  +  ykf*psif; 
iff  =  ydPpsid  +  ykPpsikd  +  yfPpsif; 

iq  =  yqq  *  psiq  +  yqkl*  psikql  +  yqk2*pstkq2; 
ikql  =  yqkl  *  psiq  +  ykql*  psikql  +y3  *  psikq2  ; 
ikq2  =  yqk2*psiq  +  y3*psikql  +  ykq2  *  psikq2; 


te  =  psid*iq  -  psiq*id; 


psidotl  =  om*psiq  -  omz*ra*id; 
psidot2  =  -om*psid  -  omz*ra*iq; 
psidot3  =  -omz*rkd*ikd; 
psidotd  =  -omz*rkql*ikql; 
psidotS  =  omz*vf  -  omz*rPiff; 
omdot  =  (omz/(2*H))*(Tm+te); 
ddot  =  om-omz; 
psidotS  =  -omz*rkq2*ikq2; 

dpsi  =  [psidotl  psidot2  psidot3  psidotd  psidotS  omdot  ddot  psidotS]'; 


SFl.M 

function  dpsi  =  sf(t,psi) 

global  ydd  ydk  ykd  ydf  yff  yqq  yqkl  yqk2  ykf  ykql  ykq2  y3  rf  rkd  rkql  rkq2  omz  ra  vf  H  Tm% 
xd  xf  xkd  xad  xaq  xq  xkq  ymd 

psid  =  psi(l); 
psiq  =  psi(2); 
psikd  =  psi(3); 
psikql  =  psi(4); 
psif  =  psi(5); 
om  =  psi(6); 
delta  =  psi(7); 
psikq2  =  psi(8); 


id  =  ydd*psid  +ydk*psikd  +  ydf*psif; 
ikd  =  ydk*  psid  +  ykd*psikd  +  ykPpsif; 


96 


iff  =  ydPpsid  +  ykPpsikd  +  yfPpsif; 

iq  =  yqq  *  psiq  +  yqkl*  psikql  +  yqk2*psikq2; 
ikql  =  yqkl  *  psiq  +  ykql*  psikql  +y3  *  psikq2  ; 
tkq2  =  yqk2*psiq  +  y3*psikql  +  ykq2  *  psikq2; 


te  =  psid*iq  -  psiq*id; 


psidotl  =  om*psiq  -  omz*ra*id  +  omz*sin(delta); 

psidot2  =  -om*psid  -  omz*ra*iq  +  omz*cos(delta); 

psidot3  =  -omz*rkd*ikd; 

psidotd  =  -omz*rkql*ikql; 

psidotS  =  omz*vf  -  omz*rPiff; 

omdot  =  (omz/(2*H))*(Tm+te); 

ddot  =  om-omz; 

psidotS  =  -omz*rkq2*ikq2; 

dpsi  =  [psidotl  psidot2  psidot3  psidotd  psidotS  omdot  ddot  psidotS]'; 


X  -10^  Fault  Current 


Time,  s 


Figure  38:  Equal  Mutuals  Model  Results  for  Krause  Simulatiou  (Phase,  Direct,  aud  Quadrature  Curreuts) 
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Figure  39:  Equal  Mutuals  Model  Results  for  Krause  Simulatiou  (Frequeucy  aud  Degrees) 
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Figure  40:  Equal  Mutuals  Model  Results  for  Krause  Simulatiou  (Field  Curreut  and  Torque) 
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APPENDIX  7  -  CALCULATION  OF  NAVSEA  DDG  PARAMETERS 


DDG-I 

DDG-II 

Phase  (KVAR) 

3125 

3125 

^ base  (^)line-line  RMS 

450 

450 

Ibase  (A)rivis 

4009 

4009 

pf 

0.8 

0.8 

Pole  Pairs 

2 

2 

H  (sec) 

2.137 

1.7083 

Ta  (Q)  p.U. 

0.00515 

0.00618 

X,  (Q)  p.U. 

0.08 

0.096 

Xaq  (fi)  p.U. 

1 

1.2 

Tkq  (fi)  p.U. 

0.0613 

0.07356 

Xkqi(Q)  p.U. 

0.329787 

0.395744 

Xad  (f2)  p.U. 

1.768 

2.1216 

rf(a)  p.U. 

0.00111 

0.001332 

Xfl(a)  p.U. 

0.136829 

0.164195 

rkd(f2)  p.U. 

0.023968 

0.028762 

Xkdl(f^)  p.U. 

0.333829 

0.400595 

Xq(a)  p.U. 

1.08 

1.296 

Xd(Q)  p.U. 

1.848 

2.2176 

Xkq(Q)  p.U. 

1.3298 

1.5957 

Xq  "(Q)  p.U. 

0.328 

0.3936 

Xd’(a)  p.U. 

0.207 

0.2484 

Xd ’’(Q)  p.U. 

0.172 

0.2064 

Xj(a)  p.U. 

1.9048 

2.2858 

Xkd(f2)  p.U. 

2.1018 

2.5222 

tdo  (sec) 

4.7846 

4.7846 

tqo”  (sec) 

0.10081 

0.10081 

tdo  ”  (sec) 

0.048521 

0.048521 

td  (sec) 

0.5553 

0.5553 

td  (sec) 

0.038912 

0.038911 
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^aq  ”1” 

Xd  =  Xad  +  Xi 


Xkq  ~  X(iq  + 


x'q=Xi  + 


Xaq  *  Xkql 
Xkql  ”1”  X(iq 


x'a^xi  + 


Xgd  *  Xfi 
Xgd  ”1”  Xfi 


Xd  = 


(  Ixad^  fxkdl^  fxf^ 


+  Xj 


Xf  =  Xad  +  Xfi 
Xkd  ~  Xdd  +  X]^dl 


1  1 


1  ,  (^r^d 

Xf 


(jOYf 


t'd  = 


V  Xgd  +XlJ  Xgd  +  XlJ 


(x)r. 


f 


(^Xkd 
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